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We study the phase transition of a real scalar <p 4 theory in the two-loop "^-derivable approximation 
using the imaginary time formalism, extending our previous (analytical) discussion of the Hartree 
approximation. We combine Fast Fourier Transform algorithms and accelerated Matsubara sums in 
order to achieve a high accuracy. Our results confirm and complete earlier ones obtained in the real 
time formalism [1] but which were less accurate due to the integration in Minkowski space and the 
discretization of the spectral density function. We also provide a complete and explicit discussion 
of the renormalization of the two-loop ^-derivable approximation at finite temperature, both in 
the symmetric and in the broken phase, which was already used in the real-time approach, but 
never published. Our main result is that the two-loop 4>-derivable approximation suffices to cure 
the problem of the Hartree approximation regarding the order of the transition: the transition is of 
the second order type, as expected on general grounds. The corresponding critical exponents are, 
however, of the mean-field type. Using a "RG-improved" version of the approximation, motivated by 
our renormalization procedure, we find that the exponents are modified. In particular, the exponent 
S, which relates the field expectation value <j> to an external field h, changes from 3 to 5, getting 
then closer to its expected value 4.789, obtained from accurate numerical estimates [2]. 
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I. INTRODUCTION 



It is well known that conventional perturbation theory fails to describe the second order phase transition of a 
tp A scalar model with Z2 symmetry in four dimensions, and more generally any system involving bosonic degrees of 
freedom, because the pcrturbativc expansion is plagued with infrared divergences [3]. This great sensitivity to the 
infrared corresponds physically to the fact that a system of bosons becomes three dimensional near a second order 
or a weakly first order phase transition. Several resummation procedures have been put forward in order to cure the 
breakdown of the perturbative expansion and thus to provide a correct description of the phase transition in a given 
model, including the order of the transition and, if the transition is of the second order, the corresponding critical 
exponents. The quality of the description of the phase transition depends usually on the level of approximation 
considered within these methods. 

The ring (daisy) resummation [4] was proposed in order to cure the infrared divergences of a massless theory through 
the resummation of the leading order thermal effects which lead to the generation of a thermal mass. This method 
produces in the effective potential a term which becomes cubic in the field at the temperature where the quadratic 
term vanishes (this is observed already at the one-loop level of the original perturbation theory) and as a result the 
phase transition turns out to be of the first-order type [5, 6]. However, it was argued in [7] that one cannot rely on 
the ring-improved perturbation theory in order to distinguish between a first- and a second-order phase transition, 
because its loop expansion parameter XT/m e f{ becomes of 0(1) at the nontrivial minimum of the potential, where 
the effective mass m e s is O(XT). 
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An even larger class of perturbative diagrams is resummed in the superdaisy [8] or self-consistent Hartrcc-Fock 
resummation scheme which, as the daisy resummation, is a local resummation scheme in that it results in a momen- 
tum independent self-energy Numerical studies and the use of the high temperature expansion revealed that this 
resummation also fails to reproduce the true nature of the phase transition [9, 10]. Recently the same conclusion was 
obtained analytically [11]. 

It was shown in [12] that in the SU(2) Higgs model the phase transition in the Higgs-Goldstone sector turns 
into second order when going beyond the super-daisy resummation by including the scalar bubble diagram in the 
propagator with which the one-loop effective potential is calculated. There are other indications in the literature that 
with the inclusion of the setting-sun diagram in the effective action the phase transition turns into a second order one 
[13-17]. Recently the type of the temperature phase transition in the ip 4 model was investigated using Monte Carlo 
simulations on a lattice [18]. It was found that for very small values of the coupling A < 10 -3 the phase transition is 
of first order, while for larger values the transition is of second order. 

The above results indicate that it is important to take into account non-daisy like diagrams in order to capture the 
nature of the phase transition. Actually, as emphasized in [19], the problems of the daisy and superdaisy resummations 
rely on the fact that, while efforts were made to include thermal effects in the effective quadratic coupling of the 
theory, nothing was done concerning the quartic coupling, which remained a constant. But in fact, as a result of 
the running, the coupling constant vanishes at T c taming around this temperature the behavior of the resummed 
perturbation theory, whose expansion parameter AT/m e ff would blow up for fixed A. An instructive comparison of 
these resummation methods with the evolution equation of the renormalization group method can be found in [20] . 
A successful description of the second order phase transition should be able to take into account the fact that the 
effective coupling constant exhibits a 4d behavior in the ultraviolet and a 3d one in the infrared. In the renormalization 
group approach the running effective coupling nicely interpolates between these two limits which allows for the correct 
description of the second order nature of the phase transition [19, 21, 22] and for the determination of the related 
critical exponents [23]. 

In this paper we study numerically the thermal phase transition of the one-component scalar field theory within the 
2PI formalism, which is known to be a systematically improvable method to resum the perturbative series [24, 25]. 
We go beyond the lowest order (Hartree) approximation used in [11] by including in the 2PI effective action the 
field dependent setting-sun diagram. This is the simplest truncation of the 2PI functional that includes nonlocal 
contributions to the gap equation for the propagator, which we solve without further approximations. In particular 
we treat the momentum dcpedcncc of the propagator self-consistently. We properly address the issue of ultraviolet 
divergences which, after being regularized using a sharp cutoff, are removed using the renormalization method recently 
developed and applied in the context of the equilibrium 2PI formalism [26-29]. 

Our main result is that the transition turns into a second order type, as compared to the Hartree approximation, 
at least for the values of the parameters that we could access. Note that another attempt to discuss the order of the 
transition from higher contributions to the 2PI effective action was pursued in [30] and in fact much more diagrams 
than the ones we shall consider here were included, namely all those which contribute at next-to-leading-order in the 
1/N expansion of the O(N) model. In this investigation however, the propagator equation which becomes momentum 
dependent at this level of approximation was solved with a momentum-independent ansatz, which lead to a slightly 
stronger first order phase than in the Hartree approximation, in disagreement with our present results which, although 
they concern the simplest nonlocal contribution to the 2PI effective effective action, involves a complete treatment of 
the corresponding momentum dependence. 

We also compute several thermodynamical quantities, as well as the critical exponents. Our conclusion concerning 
the latter is that in the two-loop ^-derivable approximations, their values remain equal to those in a mean-field 
approximation. This is not surprising since, without the inclusion of the "basketball" diagram in the 2PI effective 
action, there is no wave-function renormalization in the gap equation, and thus no possibility for an anomalous 
dimension. It is nevertheless interesting to study what happens if one implements some ideas coming from the 
renormalization group approach and let the coupling run with the temperature. We find that the values of the critical 
exponents depart from their mean field values and that some of them can be even determined analytically, such as 
the critical exponent 5 of the "magnetization" on the critical isotherm, which becomes equal to 5. As we shall see, 
this is related to the fact that, even though the approximation does not seem sufficient to generate non-analyticities 
in the field, the running coupling vanishes at the transition temperature. Some of the critical exponents have been 
studied using more elaborated truncations of the 2PI effective action by working directly in three dimensions and in 
the symmetric phase, see [31, 32]. 

Owing to the length of the text, we provide below an itemized structure of the remainder of the paper in order to 
facilitate the orientation of the reader. 

• Sec. II introduces the model, the approximation and some basic objects of the 2PI formalism which are later 
used to rcnormalize the theory. In particular, we illustrate the known fact that in a given truncation of the 
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2PI effective action, there exist two inequivalcnt expressions for the two-point function and three inequivalcnt 
expressions for the four-point function. The cutoff regularization is also introduced and discussed in details. 

• Sec. Ill motivates the need for an increased number of bare parameters when treating the renormalization of a 
given truncation of the 2PI effective action. These bare parameters are fixed by means of both renormalization 
and "consistency" conditions, which we impose at a nonzero temperature where the system is required to 
be in its symmetric phase. The consistency conditions ensure that certain features of the exact, untruncated 
2PI formalism still hold at the renormalization point and also that, despite of their increased number, all 
bare parameters are fixed in terms of only two renormalized parameters: one renormalized mass m* and one 
renormalized coupling A*. One nice feature of the two- loop approximation is that the bare parameters are given 
in terms of a finite number of perturbative sum-integrals and can thus be determined independently of the 
resolution of the gap equation. Finally some ideas borrowed from the renormalization group are implemented 
within our particular renormalization scheme. This allows us to define a RG-improved two-loop approximation 
which is later compared to the standard two-loop approximation. 

• Sec. IV presents our numerical results on the phase transition both in the two-loop and in the RG-improved 
two-loop approximation. By selecting some points in parameter space, we find numerically that where in the 
Hartrce approximation the phase transition is of the first order type, the inclusion of the setting-sun diagram 
at the level of the 2PI effective action turns the transition into a second order type, with mean-field exponents. 
The order of the phase transition is reflected also at the level of the bulk thermodynamic quantities, such as the 
heat capacity, speed of sound, and trace anomaly. We test the effects of the RG-improvement and find that the 
values of the critical exponents depart from their mean field value. 

• Sec. V is devoted to details concerning the numerical method used to solve the model in the imaginary time 
formalism. Our method is based on the fact that, in the present approximation, the self-energy only receives 
logarithmic corrections at large external momentum. As a consequence the leading part of the self-consistent 
propagator is not modified as compared to the perturbative one. Owing to this fact, to each integral involving 
the self-consistent propagator we can subtract a similar integral involving the perturbative propagator. This 
leads to a sizeable acceleration of the convergence of the Matsubara sums and to an increase in accuracy in 
the determination of convolution-type integrals with the use of fast Fourier transformations. This subtraction 
method is implemented in the gap and field equations which are solved iteratively, as well as in other quantities 
studied numerically, like the curvature and the effective potential. We investigate extensively the accuracy of our 
numerical method by testing the discretization and cutoff effects. Some related technical aspects arc discussed 
in App. B, where a collection of perturbative integrals is also given, and in App. C, which is devoted to the 
acceleration of the Matsubara sums. 

• Sec. VI shows that the expressions for the bare parameters obtained from the renormalization and consistency 
conditions render finite the gap and field equations, as well as the effective potential. Some more technical 
parts are relegated to App. A. This section has only theoretical relevance in that it shows what needs to be 
done in order to check explicitly the renormalizability of the model in the present approximation, but the finite 
equations derived here are not used to solve the model numerically. 

• Sec. VII is devoted to some conclusions and outlook. 



II. THE TWO-LOOP ^-DERIVABLE APPROXIMATION 

In this paper, we consider a real scalar tp theory in four dimensions at finite temperature, defined by the Euclidean 
action 

%] dr J d*x Q(cW) 2 + i(V^) 2 + !=jv + , (1) 

where the inverse temperature sets the range of integration over the imaginary time. The parameters m and Ao 
denote respectively the bare mass and the bare coupling. To ensure that the spectrum of the underlying Hamiltonian 
is bounded from below, one should restrict in principle to positive values of Ao- We shall discuss this condition in 
more detail as we treat renormalization in Sec. III. 
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FIG. 1: Two-particle irreducible diagrams contributing to <&[</>, G] in the two-loop ^-derivable approximation. Plain lines 
represent the propagator G and circled crosses represent the field 0. 



A. Effective potential and gap equation 

The two-particle-irrcduciblc (2PI) formalism provides a representation of the effective potential 7(0) in terms of 
2PI diagrams. It is obtained as the value taken by the functional 

7 [0, G] = -^-<t> 2 + "jj"0 4 + \ ( T [ In G- l (Q) + (Q 2 + m 2 )G(Q) - l] + <I>[0, G] (2) 

J Q 

at its stationary point G = G, that is 7(1/)) is equal to 7 [0,(5] with = <57[0, GJ/JGIg. In Eq. (2), the variable 
represents a homogeneous field configuration and G(Q) = G(iuj n ,q) an even function of the Matsubara frequency 
u> n = 2irnT and of the modulus of the three dimensional momentum q = |q|. We have also adopted the notation 

/ f(Q) = TJ2 f(i"n,q) = TJ2 hAsfi^q)- (3) 

J Q n J 1 n=-oo J W 

The functional derivative with respect to G, which appears in the definition of G, is to be understood as the functional 
derivative in the space of functions G(iu n , q) which are even with respect to u n and rotation invariant with respect to 
q. However, in all expressions where such functional derivatives appear in this paper, they can be treated as normal, 
that is unconstrained, functional derivatives. Note also that our definition of the functional derivative in Fourier 
space implies a factor (2ir) 3 /T. With this convention, if J-[G] is a functional of G which we evaluate for G = G, the 
following chain rule applies 

clF[G] _ f T SF[G] dG(Q) 



90 J Q SG(Q) 



G 



Finally, the functional —$[0, G] corresponds to the sum of all 0-leg 2PI diagrams that one can draw in the "shifted" 
theory S[<j> + (p] — S*[0] — (5S/5(j))ip at finite temperature, using the function G in place of the free propagator. This 
functional cannot be computed exactly So-called ^-derivable approximations consist in retaining in $ [0, G] only 
certain classes of diagrams. In this paper, we consider the two-loop ^-derivable approximation: 

$[0,G] = ^0 2 T[G] + ^T 2 [G] - |0 2 5[G] , (5) 
which corresponds to the 2PI diagrams represented in Fig. 1. We have introduced the notations 

T[G] = f G(Q) and S[G] = f f G{Q)G(K)G(K + Q) (6) 

JQ JQ JK 

for the "tadpole" and "setting-sun" sum-integrals, to be used throughout this work. Similarly, we introduce the 
notation 

B[G]{K)= ( T G(Q)G{Q + K) (7) 



for the "bubble" sum-integral. The setting-sun sum- integral reads then S[G] = Jq G(Q)B[G](Q). 

According to the above discussion, in order to evaluate the effective potential, one needs first to determine the 
propagator G from the stationarity condition = 5^[4>,G]/SG\q. This condition can be rewritten as G~7 1 T {Q) = 
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M IAQ) with Q 2 = + 1 2 and 



(8) 



which we refer to as the "gap equation" . We have used momentarily the subscripts <f> and T to stress the fact that the 
propagator G^ ; t(Q) and the momentum dependent mass M^^iQ) depend both on the field </> and on the temperature 
T . In what follows, we shall omit this notation unless specifically needed. A particular role will be played by the 
value of M 2 (Q) at Q = 0, which we denote more simply by M 2 . In the two- loop ^-derivable approximation, the gap 
equation reads 

M 2 (K) =m 2 + ^ + ^ T[G] - f 2 B[G](K) , (9) 
which we obtained from Eqs. (5) and (8) by making use of the identities 

= > - i§ = ™' (10) 



B. Field equation and geometry of the effective potential 



The phase transition will be studied by monitoring the position of the extrema of the effective potential as the 
temperature T is lowered from an initial temperature at which the system is required to be in its symmetric 
phase, see the discussion in Sec. IV, down to zero temperature. The "field equation", which codes the position of 
the extrema, is easily obtained if one makes use of the stationarity condition = Sj[4>, G]/5G\q to express the first 
derivative of the effective potential as 



5 1 5 7 [0,G] 



m o 0- 



6 ' 



d$[0,G] 



A 

V 



Ao 



G 



2 T[G]-^S[G] 



(11) 



Note that the field equation, which is obtained by equating this first derivative with zero, is coupled to the gap 
equation. Thus, for the purpose of determining the extrema of the effective potential, the gap and field equations 
need to be solved simultaneously. 

Some valuable information can also be obtained from the field derivatives of the effective potential, in particular 
from the second and fourth derivatives at <p = 0, at least while the potential is defined around = 0. Taking a second 
derivative with respect to in Eq. (11) and evaluating the result for — 0, we obtain 



6^2 

A^2 



A, 



\2 



S[G^ =0 ] 



(12) 



0=0 



We shall later use this formula in order to define and determine the critical temperature in those cases where 
the system undergoes a second order phase transition. Note that we could obtain an expression for the second 
derivative M 2 = 5 2 j/S(j) 2 for any value of the field but it is substantially more complicated than Eq. (12). Moreover, 
the definitions of M 2 and M 2 generalize in a straightforward way to any ^-derivable approximation and, if no 
approximation was considered, these two quantities would coincide. The fact that M 2 and M 2 are not equal in a 
given approximation needs to be regarded as a truncation artifact: here for instance, M|_ — M 2 =0 = C(A§) is a 
discrepancy that lies beyond the accuracy of the present approximation. This fact, and a similar one concerning the 
four-point function that we discuss now, is the main motivation for the renormalization procedure that we present in 
Sec. III. 



Similarly, we can take three field derivatives on Eq. (11). Evaluating the result for = 0, we obtain 



<5 4 7 



A + 3 



0=0 



f 


"a 5T[G) 


X 2 5S[G] 




d 2 G 2 (Q) 




Iq 


2 SG{Q) 


Co 6 6G (Q) 




dej) 2 


4>=o 



(13) 
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where we have used Leibniz rule for computing multiple derivatives of the product of two functions, hence the factor 
of 3, and the chain rule (4) together with dG/d<j>\^ Q = 0. From Eq. (10) and d 2 G / 'd<j)%= = -G\^d 2 M 2 / 'd(t)% =0 , 
we obtain finally 



v ^ 



= 



A - Aoi3[G0 =o ](Q) G 2 =0 (Q) — ® 



(14) 



= 



The quantity d 2 M 2 /d4> 2 \^, = o obeys a linear integral equation which can be obtained from the gap equation (9) using 
the same strategy as the one that leads to Eq. (14). We obtain 



d 2 M 2 {K) 



X -X 2 B[G^ =0 ](K) 



Ao 



1=0 



G 2 (O) d ^ PiQ) 



1=0 



It will be convenient to introduce some additional notations. We shall write Eq. (14) as 



A0=o(Q)G 2 =o (Q)^ =o (Q), 



with A^ = q = Ao, A^=o(K) = A - A 2 , B[G^=o](K) and 



d 2 M 2 {K) 



(15) 



(16) 



(17) 



0=0 



The quantity V^ = o(K), whose value at K = we denote more simply by V^—o, obeys the linear integral equation (15) 
which we rewrite as 



V^K) = K^ =0 {K) 



A0=o 



with A^, = o = Aq. 1 This equation can be solved explicitly as 2 



in terms of the quantity V^ = o such that 



^0=0 



Gl =0 {Q)V^ =0 {Q) . 



G| =o (g)A 0=o (O), 



(18) 



(19) 



V*= 



= 



Afc 



= 



S[G0= O ](O) 



(20) 



As it was the case for M^ =0 and M 2 =0 , the definitions V^=o, V^>=o and V^=o can be generalized to non- vanishing values 
of the field and can be defined for any ^-derivable approximation. Moreover, when no approximation is considered, 
they coincide with each other and represent a unique definition for the four-point function at zero external momentum. 
This is one of the clues to understand the renormalization of ^-derivable approximations [28] which we illustrate in 
Sees. Ill in the case of the present two-loop approximation at finite temperature. 



1 The reason for introducing two different notations A^,— o and A^,— q shall become clear in Sec. III. 

2 To see this, one writes first Eq. (18) as 



f 

■!() 



A, 



=0 /=,2 



where (K — Q) is the product of a Kronecker delta <5 n m for the frequencies, and a three dimensional Dirac distribution <5( 3 '(k — q). 
This equation can then be easily inverted in terms of V^—q because the definition (20) is equivalent to 



S W(P_ K) _^ (if) 



S^(P-Q). 
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C. A few words on regularization 



The equations derived in the previous sections will be used throughout this work but strictly speaking they do not 
make sense in the absence of an ultraviolet regularization. Before any practical application, it is therefore mandatory to 
give them a precise meaning by choosing some regularization. To obtain a proper regularization of the 2PI functional, 
we start from the functional W[J, K] defined by 



T>ip exp 



1 



<p ■ (Go^a)" 1 -<p > x 



/ Dip exp • {G RaY 



S] r 



J- 



9- 



K-<p} 



JVtpexp{-^tp- (Gq-Ra)- 1 -<p} 



(21) 



f = f x J ( x )<P( x ) and <P - K • <p = f x J <p(x)K(x,y)<p(y), with J x 



where we have introduced the notations J 

Jq^ T dr J d 3 x. Wc have written the quadratic part of the action ip ■ (Go-Ra) -1 ■ f in terms of a regulating function 
R\ which in 3d-momentum space cuts off momenta q > A and obeys the property Ra{<1 < A) -> 1. In this work we 
shall restrict to a sharp regulating function R\(q) = 0(A — q). However, in the remainder of this section, we consider 
an arbitrary regulating function. We have also introduced a normalization such that the second factor of Eq. (21) is 
completely regularized by the presence of Ra , at least perturbatively and for sources K close to K — (this can be 
checked by expanding the second factor perturbatively since all the Fcynman diagrams involve the fastly decreasing 
propagator GqRa)- The first factor of Eq. (21) requires its own regularization but this is straightforward since this 
factor is Gaussian. It can be written cxp(— 70 (too, A)) where 7o(too, A) corresponds to the free energy per unit volume 
of the free theory (obtained for SmtM = 0): 



7o(m , A) 



d 3 q 
(2tt): 



Ra(q) 



,(0) 



2Tln(l 



'An 



(22) 



with £q = \J q 2 + TOq. After Legendre transformation of W[J, K) and restriction to homogeneous field configurations, 
one obtains the "regularized" functional 



7 [^G]= 7 o(TO ,A) + ^^ + ^ 4 + i 



lnG _1 (Q)-La- 



Q 2 



R A (q) R A (q) 



G(Q)-1 



+ mG}. (23) 



What is meant by a regularized functional is not completely straightforward. In fact, a given functional of G can be 
well defined for certain classes of propagators but not defined for others. The important fact about the functional 
(23) is that it is well defined for a class of propagators in the vicinity of its stationary point. 3 In order to check this 
statement, note first that the gap equation, which defines the stationary point, reads 



G~\Q) 



Q' 2 



Ra(q) 



n(Q) with n(Q) 



2(5$ 



SG{Q) 



(24) 



One can convince oneself that, due to the presence of Ra(q) and even though the Matsubara sums are not cut off by 
the regularization, any nonlocal contribution to the self-energy Tl(Q) is suppressed 4 at large Q. Then, in this limit, 
the self-energy n(Q) approaches a constant ll^ which is entirely determined by the first two diagrams contributing 
to $[</>, G] (the only ones which give a local contribution to the self-energy), see Fig. 1: 



fl(Q) -> floo = 



A 



Ao 



T[G] 



\Q\ 



(25) 



It follows that, at large Q: 



G(Q) 



Q 2 + m 2 (Q 2 + m 2 ) 2 



(26) 



3 More rigorously stated, the functional has an extremum within the class of propagators for which it is defined. This is in one to one 
correspondence with the fact that the functional (21) is regularized in the vicinity of K = 0. 

4 For instance, if one considers the perturbative bubble diagram B[G](K), after performing the Matsubara sum in Eq. (B15) of App. B, 
it is clear that the regularized diagram behaves like 1/lu 2 at large uj and vanishes if |k| > 2A. 
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where the second term is sublcading with respect to the first one. From this we deduce first that all the diagrams 
appearing in the gap equation or in the effective potential through the functional <&[</>, G] are regularized, which was 
not obvious a priori. Moreover, the explicit sum-integral appearing in the 2PI functional (23) is well defined when 
evaluated for G = G. This is because at large Q we have: 

For propagators G in the vicinity of G, we expect this behavior to be only slightly modified, which shows that the 
regularized 2PI functional (23) is well defined in the vicinity of G = G, as announced above. 

For practical purposes, it is convenient to consider the change of variables 5 G — > GRa- In terms of this new variable, 
the regularized 2PI functional becomes 

7 [0, G] = 7o (mo, A) + + ^ + l - J T R A (q) [ In G~\Q) - ln(Q 2 + m 2 ) + (Q 2 + m 2 )G(Q) - l] + <f>[0, GR] , 

(28) 

where we have introduced an additional Ra{q) in the explicit sum- integral of Eq. (28), although it is not needed 
a priori due to Eq. (27). This is convenient however because the gap equation takes then the simple form (for 
i?A(g) = 0(A — q), this equation needs to be understood for q < A) 



G-\Q) = Q 2 +rrr 



8G(Q) 



(29) 



GR 



Moreover, the first derivative of the effective potential 7(0) = 7 [0, G] reads 

(30) 



GR 



It follows that the formal expressions obtained in the previous sections for the gap and field equations and also the 
different n-point functions can be regularized by replacing each propagator G by GR. From now on, we will thus 
consider that such a replacement has been done, but we shall leave the regulating function implicit. As far as the 
effective potential is concerned, it will be computed by evaluating the functional (28) for G = G. In fact, it is more 
convenient to use the following formula (m 4 denotes for the moment an arbitrary parameter) 



(31) 



to rewrite the regularized effective potential (28) as 

7 [</>, G] = 70 (m*, A) + ^ +^ + \ f R ^V* G ~\Q) - HQ 2 + ml) + (Q 2 + m 2 )G(Q) - l] + $[0, GR] . 

(32) 

As a final remark, note that the rcgularization that we have introduced, because it only cuts the modulus of the three- 
dimensional momentum, breaks explicitly the four dimensional rotation symmetry of the theory at zero temperature. 
However, since the operators (d T ip) 2 and (V</?) 2 do not require any renormalization in the present approximation, see 
below, the continuum results obtained as the cutoff is sent to infinity possess the four dimensional rotation symmetry 
at zero temperature. 



5 Strictly speaking, if we want the exact identities M 2 = M 2 and V = V = V to hold true after this change of variables, we should also 
transform the field according to <j> — > (py/RX in momentum space. Note however that since we restrict our study to homogeneous field 
configurations, the field in momentum space is concentrated around its zero mode component which is not affected by the change of 
variables since by assumption the regulator is such that Ra(<? < A) -> 1. 

6 At least as long as we restrict to homogeneous field configurations. 
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III. RENORMALIZATION PROCEDURE 

As in the case of the Hartree approximation [11], we will show that it is possible to adjust the dependence of the 
bare parameters in such a way that the results become insensitive 7 to the regulating scale A as the latter is taken to 
infinity. In contrast to perturbation theory, and due to certain truncation artifacts which appear within ^-derivable 
approximations, we will need to introduce more bare parameters than usual. An important step will then be to 
fix all these parameters in terms of the usual number of renormalized parameters: one renormalized mass and one 
renormalizcd coupling. 

A. Multiply denned bare parameters 

A general remark is in order first. The two-loop approximation is not rcnormalizable a priori in the form we 
have presented it so far. To see why this is so, let us consider for instance the quantities M%_ and M%_ Q . We 
have already seen that these two quantities differ by contributions of order Aq. It follows that their quadratic 
divergences differ by terms of the same order. On the other hand, a closer look at Eq. (9) for <f> = and 
Eq. (12) shows that the bare mass squared ttiq enters both equations in exactly the same way, that is as a tree 
level term. Then, M| =0 and M| =0 cannot be renormalized simultaneously using the same bare mass: if one 
adjusts Trip to absorb the quadratic divergence in one of the masses, there is an unbalanced quadratic divergence 
in the other one and vice-versa. Similar remarks apply to the different definitions of the four-point function 
at zero external momentum, V^q, V^ = o and which, although they differ at order Ag, involve the bare 

coupling Ao at tree level in an identical manner, leading once more to unbalanced divergences. Two attitudes are 
possible from this point on: one can either consider that ^-derivable approximations are ill-defined regarding the 
issue of renormalization, or one can try to cure these truncation artifacts and define a sensible renormalization scheme. 

There is actually a simple solution to the problem of unbalanced divergences: in the case of the masses, one can 
slightly modify one of the two equations, Eq. (12) for instance, by introducing a second bare mass parameter TO2 
in place of mo- Similarly, one can replace the tree level contribution Ao by A2 in V^ = o and by A4 in V^ = o- Despite 
the apparent simplicity of this way out, it is important to bear in mind that this procedure is only acceptable if it 
provides a way to determine ni2, A 2 and A4 without introducing more physical or renormalized parameters than the 
ones which are usually present in a scalar tp 4 theory and if it ensures that the discrepancy between mo and m-2 and 
between Ao, A2 and A4 disappears as the order of truncation is increased. We shall treat these matters in the next 
subsection. For the moment, we just note that all these changes can be formulated at the level of the regularized 2PI 
functional (32) which reads now 

7 [0,G] = 7o ( m „A) + ^ 2 + ^0 4 + i^ [\nG- 1 {Q)-HQ 2 + ml) + {Q 2 +ml)G{Q)-l] + *[<t>,G\, (33) 
with 

m G] = ^ 2 T[G] + ^T 2 [G] - ^ <j?S[G] , (34) 

and where, as explained in the previous section, each propagator G(Q) in each diagram contributing to G] 
is multiplied by R\(q) = 8(A — q) and the explicit integral over q in Eq. (33) is cut off at the scale A. Note 
that there is an additional modification which we have not mentioned so far: the bare vertex which appears in 
the highest loop diagram, that is the setting-sun diagram, has been replaced by A*, which will be identified below 
as the renormalized coupling at some renormalization scale. This replacement is again necessary if one wants 
to avoid unbalanced divergences: keeping the bare vertex in the setting-sun diagram would absorb divergences 
related to 2PI diagrams which are not present in the two-loop approximation. Actually, this is also what occurs 
in a perturbativc calculation. However, in contrast with the perturbative case where replacing Aq by A* in the 



7 As explained in [33] , the discussion of divergences in a non-perturbative context such as the one considered here might be different from 
that in perturbation theory. In particular, certain divergences of the perturbative expansion do not appear as such after resummation. 
However they lead in general to the same difficulty, namely to the fact that one cannot define cutoff insensitive results with a high 
accuracy unless renormalization of the mass and the coupling is considered. Throughout this work the term "divergence" will be used 
in this somewhat extended meaning. 
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setting-sun diagram can be interpreted in terms of an expansion in powers of A*, the situation is more subtle 
in the case of ^-derivable approximations and one needs to provide a consistent and systematic way to fix the 
bare coupling which appears in the setting-sun diagram and more generally in higher loop 2PI diagrams. Such a 
procedure exists and is pretty similar to the one we shall present below for fixing the additional bare couplings A2 
and A4 but it lies slightly beyond the scope of the present paper. In what follows, we shall then admit that the re- 
placement Ao — > A* at the level of the setting-sun diagram is consistent at this order and refer to [28] for further details. 

With all these modifications taken into account, the gap equation becomes 

M\K) = mg + ^ + £ T [G] - I ^ B[G]{K) , (35) 
and the first derivative of the potential reads 

g = ,(^ + ^ + | riei _| s|ei ), m 

from which one obtains the field equation 

= (mj + y 2 + y T[G^\ - 1 . (37) 
The expression for the second derivative of the effective potential (the curvature) at </> = is 

M| =0 = ml + ^ T[G^ ] - ^ S[G^ ] . (38) 

Finally, the different definitions of the four-point function, that is Eqs. (20), (18) and (16), remain unchanged if we 
replace the previous definitions of A^o, K^, = q{K) and A^o by 8 

A 0=o = Ao, K^{K) = \ 2 -\lB[G^]{K), and A 0=o = A 4 . (39) 

Our results regarding renormalization concern this precise formulation of the two-loop <I>-derivable approximation. 
We will show that it is possible to fix the dependence of the bare parameters too, 771.2, Ao, A2, and A4 with respect to 
the scale A such that the solutions to the gap and field equations, and the effective potential converge as A is sent to 
infinity. The proof is quite technical, in particular due to the fact that one needs to prove that the bare parameters 
can be taken independent of the field cj> and of the temperature T. For this reason, we postpone the proof until 
Sec. VI, where the interested reader can find all the details. For practical purposes, we only need to know that such 
a proof exists. The expressions for the bare parameters can then be obtained by imposing appropriate conditions, as 
we explain in the next two subsections. Note however that all the conditions cannot be independent because, despite 
the higher number of bare parameters required to absorb the ultraviolet divergences, the theory is to be parametrized 
in terms of the usual number of physical or renormalized parameters. 

B. Parametrization in terms of two renormalized parameters 

Allowing for more bare parameters than usual is only acceptable if these are fixed in terms of the usual number of 
renormalized parameters and also if the renormalization procedure is such that the discrepancies between TO2 and too 
as well as between A4 , A2 and Ao are pushed to higher orders as one increases the number of loops of the ^-derivable 
approximation. In order to understand how this is possible within a given ^-derivable approximation, the crucial 
point is to remember that the need for multiple bare parameters arises from the existence of multiple definitions for 
the n-point functions, which coincide when no approximation is considered at all. We have seen in particular that 
it is possible to define the inverse two-point function at zero momentum in two different ways and the four-point 
function at zero momentum in three different ways. Now, if we imagine for a moment that these quantities were 
measurable at some temperature T*, there would be only one measured value for the mass at this temperature, and 



Note that A^,— q and A^—q become two different quantities, which explains why we introduced two different notations in the first place. 
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only one measured value for the coupling at this same temperature, in apparent contradiction with the multiplicity 
of definitions for the two- and four-point functions. It is then quite natural to adjust the different bare parameters 
in such a way that these truncation artifacts disappear at the temperature where we would make contact with a 
measurement. 

All the bare parameters that we introduced in the previous section can then be fixed through the conditions 

M 4,=o, t, = M l=a, 7\ = m l and %=o, T t = V^o. Ti , = %=o, t* = A* . (40) 
Note that we can arrange the previous conditions in two categories. Two renormalization conditions, for instance 

M*=o,t. = m l and V*=o,T, (41) 

and three "consistency" conditions 

M l=o,T* = M l=o,T„ , ^=o,T, = Vfi=o, Tt , and %=o, Tt = Vfi=o, Tt ■ (42) 

The consistency conditions do not involve any renormalized parameter. In this way, it is possible to fix the bare 
parameters in terms of the usual number of renormalized parameters. These conditions ensure also that the discrep- 
ancies between mo and mi or between Ao, A2, and A4 are beyond the accuracy of the approximation at hand, as we 
will check on the explicit expressions for the bare parameters that we obtain below. For practical purposes, we shall 
use the consistency conditions in the form 9 

M 4>=o, t, = m l > v <t>=o, t„ = A* , and %=o, t, = A* , (43) 

which arc obtained trivially by combining Eqs. (41) and (42). Note finally that m 2 needs to be taken positive in 
order for the gap equation at (ft = and T = to make sense. From the consistency conditions, it follows then 
that Mj_o t„ * s positive: the effective potential is thus convex around (ft = at the renormalization temperature 
T*. We will see later that the effective potential is in fact globally convex at this temperature. Our system is then 
parametrized in the symmetric phase. We shall also see below that the sign of A* needs to be taken positive. 

C. Explicit expressions for the bare parameters 

We now use the renormalization and consistency conditions in order to determine the expressions for the bare 
parameters mo, 771 2, Ao, A2, and A4. As we show in Sec. VI, these expressions are such that the solutions to the gap 
and field equations, as well as the effective potential (up to a field and temperature independent constant) become 
insensitive to A at large A. 

From Eq (35) and the condition (41) for M 2 , we obtain 

ml = m% + ^-T4G*}, (44) 

where 7i[G*] = Jq* G?*(Q+) and G*(Q*) = l/(Ql + M^ =0 T J = l/(Ql + ml) is the propagator at temperature T* and 
eft = 0. The notation is used to emphasize the fact that the Matsubara frequencies are considered at the reference 
temperature that is = (iw*, q) with u>* = 2imT+. Since G*(Q 7t ) is a free-type propagator, we obtain an almost 
explicit expression for mo in terms of m*: 

m§ = «£-yT*[G*]. (45) 



Another interesting possibility is to impose the consistency conditions (42) at a temperature T* different from T*, while using the 
same renormalization conditions (41). This introduces TVdependencies which in some sense can be interpreted as scheme dependence. 
They could then be used to test the convergence of the <I>-dcrivable expansion by comparing the sensitivity of the results to the scale 
T* between two orders of approximation, for instance between the present two-loop approximation and the three-loop approximation. 
Although interesting, this is beyond the scope of the present work. 
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From Eq. (38) the condition (43) for M 2 , we obtain similarly 



m% = ml - yT*[G*] + ^S*[G*] , 



(46) 



where <S*[G*] = j£ G^Q^G^K^G^K^ + Q*). We verify below that A 2 - A = 0(A 2 ) from which it follows 



that 



0(A 2 ), that is the discrepancy between the two bare masses is beyond the order of the present 



approximation. This feature generalizes to higher order truncations and is related to the use of the consistency 
conditions. 10 



We proceed similarly for the bare couplings. The condition (41) for V^ = o leads to 



(47) 



where B*[G*](0) = /^* G 2 (Q*). As the cutoff A increases, there is a scale at which Ao diverges and above which it 
becomes negative. This is the Landau scale A p defined by the condition 



° = a; 



e(A p - q) 

Qt {Ql + mlY 



(48) 



where we have made the bubble sum- integral and the ultraviolet regulator explicit. If one wants to maintain Ao 
positive, one needs to choose A* positive and A below the Landau scale A p . We shall work with values of the 
parameters such that A p is much larger than all other scales in the problem, namely to* , T* , T, <fi <C A p . Then our 
results will be pretty much insensitive to the cutoff A in the regime to*, T*, T, <f> <C A < A p . Note that, in the present 
two-loop approximation, it is mathematically possible to take A > A p at the level of the renormalized quantities, 
and even consider their continuum limit as A — > oo. The difference between this continuum values and the values 
obtained for A < A p are pretty tiny. 11 

The bare coupling A2 is determined from the condition (43) for V^—o and from the condition (41) for V^—o- Using 
Eqs. (18) and (39) we obtain 



A* = A 2 -A 2 ^[GJ(0)-^ / Gl(Q*)[\ 2 -\lB4G+}(Q+)] 



It is convenient to decompose A2 as A 2 = A21 + 6\2ni with 

SX 2 nl = AX[G*](0) . 



(49) 



(50) 



It follows that 



A21 
A21 



1 



A, 

A2 



Ao A * + 2 



S*[GJ(0) 
Q 



A| 
2 



G;(Q*)[B*[G*](Q*)-B*[G*](0)] 



G 2 (O,)[^[G.](Q,)-^[G,](0)], 
where we have used Eq. (47). We arrive then at 

1- y / *G 2 (Q*)[S*[G*](Q*)-R[G,](0)] 



(51) 



(52) 



In fact there exists a systematically improvable class of <t>-derivable approximations for which M^_ T = M^_ T and thus m| = rrip, 

sec the discussion in [28]. Within this class of truncations, one has also V^,— o,t = V<t>=o,T an( i thus A2 = Ao. However V^ = o t V#=o,T 
and thus A4 ^ Ao- 

In approximations where a Landau pole could appear as an actual pole of the propagator, the continuum limit does not exist. However, 
the effect of renormalization can still be understood as an insensitivity with respect to the cutoff, up to terms of the order of inverse 
powers of A, in the regime m*, T* , T, <f> <C A A p . 
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Finally, the condition (43) for V^—q given in Eq. (16) with A^ = A4 leads to 



A 4 = A* + - 



V*(Q*)G*(Q*)A*(Q*). 



(53) 



where 12 



Then 



A*(iQ = A^o.tAK*) = A21 - A2[B*[G*](if*) - B*[G*](0)] , 
F*(#*) = ^= ,T.(«;) = A*-A2[iB*[GJ(Jif*)-B*[G!*](0)]. 



A 4 = A* + ^A*A 2 i£*[Gy(0)-^(A 21 + A*)A* / * G 2 (Q*) [<B*[Gy(Q*) - B*[G*](0)] 
+ ^A 4 { T * Gl(Q*)[B*[GA{Q*) - B4G*}(0)] 2 . 

2 .In 



(54) 
(55) 



(56) 



Using Eqs. (47) and (52), this becomes 



and thus 



A 4 = A, + 3-^(A - A*) + 3(A 21 + A,) - 1 
Ao V *o 



3 r T* 

A 4 
2 * , 



G£(QO[B*[G*](QO-B*[G*](0)] 2 , 



(57) 



_21 

A 



G*(QO[B*[G*](QO-B*[G*](0)f 



(58) 



Note that X 2 — A = O(Aj), as announced above. Similarly A 4 — A = O(A^). Moreover, since Ao increases strictly 
with A in the interval [0, A p ] from the value A* at A = + to 00 at A = A+, we deduce using Eqs. (50), (52) and (58) 
that the bare couplings A 2 and A 4 are positive for A < A p , that they are equal to A* at A = + and that they diverge 
for A = A+. In the case of A2, this uses the fact that B*[G*](Q*) — Z?*[G*](0) is negative, which we prove in App. B. 
Thus, in the present approximation, although we had to introduce different bare couplings to cope with the problem 
of unbalanced divergences, the Landau scale at which all these bare couplings diverge is uniquely defined. 

All the bare parameters involve perturbative sum- integrals in terms of the free-type propagator G* . In fact all the 
bare parameters but A 4 can be reduced to the tadpole, bubble and setting-sun perturbative sum-integrals (and their 
mass derivatives), which makes their evaluation relatively easy: the Matsubara sums are computed exactly and the 
remaining integral over the modulus of the momentum is determined using adaptive integration routines. This is 
obvious for mo, 7712 and Ao and <5A2 n i- For A21, we can use the formula 



A, 



(59) 



In the case of A 4 , we can reduce it to 



A 4 = A* + 3-^ - 3A* I 3 - 2^ 



Ao 



Ao 



-a* 

2 * 



r. 



Gl{Q*)Bl[GA(Q*) 



(60) 



The last term cannot be reduced to simpler sum-integrals and is then evaluated directly as a double sum. The 
variation of the different bare couplings with the cutoff A is shown in Fig. 2, for two different values of the 
renormalized coupling A*. In the left panel, corresponding to a large value of the coupling A* = 8, the Landau scale 
is relatively close to the scales T* or m*. We observe the divergence of the three bare couplings at the Landau scale. 



We use the fact that V^ =0 (K) - V</, =0 (0) = K^ =Q {K) - A^ =0 (0). 
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FIG. 2: Variation of the three bare couplings Ao, A2, and A4 with the cutoff A, for m*/T* = 0.04 and T* = 1. The left panel 
corresponds to a renormalized coupling A* = 8 for which the Landau scale is A p /T* ~ 540. The right panel corresponds to a 
renormalized coupling A* = 3 for which the Landau scale is A p /T* ~ e 39 . The lines are obtained by performing exactly the 
Matsubara sum and evaluating the integrals over the modulus of the momentum using adaptive integration routines, except 
for the last term of Eq. (58) which is evaluated as a double sum. The points are obtained by evaluating the integrals in the 
expressions of the bare couplings as a double sum using N T = 2 10 non-negative Matsubara frequencies and N s = 3 x 2 10 values 
of the modulus of the the 3d-momentum. The discrepancy between the points and the corresponding line is related to the 
discretization of the momentum integrals. This will be used in Sec. V in order to discuss discretization effects. Note also that 
the Matsubara sum in the expression of Ao was accelerated using Eq. (C12). 



In the right panel, corresponding to a small value of the coupling A* = 3, the Landau scale is pretty far apart from 
the scales T* and m*. For cutoff scales A in the regime to*,T* « A < A p , we observe a logarithmic dependence of 
the bare couplings which mirrors the divergences of the gap and field equations absorbed by these bare parameters. 
In particular, the variation of A4 is rather important and shows that, if not properly renormalized, logarithmic 
divergences can lead to sizable cutoff dependencies, even if the coupling is small. 

As a final remark note that starting from Eq. (33) one can derive an expression for the difference Aj(<f>) = 
7(0) ~~ 7*(0), where 7*(0) is the value of the potential at temperature T* and at vanishing field. 13 Hereinafter this 
difference will be referred to as the subtracted effective potential. To derive it, we use the expression (36) for 8j/5<fi 
as well as the expression (45) for m§. Noting that 7*(0) = 70 (m*, A) — (An/8)7^ 2 [G*], it is then simple to arrive at 
the following expression: 

In G~ X {Q) - lnG-^Q) + (Q 2 + m 2 )G(Q) - 1 

- ^ + y^ + ^ r[Q] - %[G ^ 2 - (6i) 

As far as the nature of the transition is concerned, it is enough to concentrate on A'j^cj)). The expression above is 
useful on a practical level because /5(f) needs to be computed anyway when solving the coupled system of gap and 
field equations and we can then use the same numerical routine. Moreover, Eq. (61) involves a difference of tadpole 
sum-integrals which can be computed efficiently, as we explain in Sec. V. We will also prove in Sec. VI that Aj^cj)) is 
finite, whereas 7(0) is only finite up to a temperature and field independent divergent constant. 



A7(0)=7o(m*,A)-7j(m*,A) + 



D. Renormalization group improvement 

In addition to solving the two-loop ^-derivable approximation, we shall also consider an "improved" two-loop 
approximation based on some ideas borrowed from the renormalization group and which we now explain. 



This quantity is well defined because according to the renormalization condition for M^_ , M^_ T is defined for T = T*. 
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In the renormalization procedure that we have presented in the previous sections, the temperature T + played 
the role of a renormalization scale [i. In the exact theory, the physical observables should not depend on /i: any 
change in fj, should be compensated by a "running" or "flow" of the renormalized parameters m*(/x) and A*(//). In 
principle, if one is able to determine the running of the parameters, it is then possible to describe the same theory 
from different but equivalent points of view, each implying its own renormalization scale and the corresponding 
renormalized parameters. In particular, in calculations at finite temperature, one can choose a description in which 
the renormalization scale equals the temperature T. 

The previous considerations become particularly interesting in the presence of some approximation because the 
different possible descriptions cease to be strictly equivalent. It can then happen that taking into account the running 
of the parameters leads to an "improved" approximation. Usually, the improvement is related to the fact that the 
running resums higher order contributions. In the present work, we shall see that the running will have somehow the 
opposite effect in the sense that it will remove certain fluctuations, namely fluctuations responsible for some of the 
artifacts of the ^-derivable approximation that we mentioned above. 



In order to obtain the running of the renormalized parameters with the scale T in the present approximation, we 
choose Eqs. (45) and (47) and differentiate them with respect to T+ under the assumption that the bare parameters 
mo and Ao are fixed. 14 Then, m*(T) and A*(T) can be obtained by integrating the ordinary differential equations for 
dA*(T*)/dT+ and dm1(T+)/dT+, starting from the initial temperature T* at which we fix the value of the renormalized 
parameters: m^(T+) = and A*(T+) = A*. In the present approximation, there is in fact an easier way to proceed. 
Indeed, by comparing Eqs. (45) and (47) with 

M$=o = < + Y^=o], (62) 
Js[G 0=o ](O), (63) 



Vrf, = o A 2 

we see that, since M0 =o ,t 4 = m * = m *{T±) and V^o.t* = A* = A^(T^), the dependence of m*(T) and A*(T) on T is 
nothing but that of M^ = o and V0=o on T. This simple fact provides us with the following recipe to implement the 
RG-improvement : 

1. solve the gap-equation (62) for Af^ =0 in terms of the parameters T+, m+ and A*; 

2. compute V^o from Eq. (63), using the determined M^, = o; 

3. apply the replacements — > T, —¥ M^—o, A* — } V^—o in every equation of interest. 

The replacements apply also to the bare parameters m 2 , A 2 , and A4, which have to be redetermined and will be 
denoted m^, A^* 3 , and Af G when needed. The bare parameters nig and Ao do not need to be modified since they are 
invariant, by construction. 

As an illustration of how the improvement works, let us consider the curvature of the effective potential. Before 
the improvement, it reads 

M|=o = m l + Y [ T ^=°] - T *[ G *]] - y [SiG^o] - <5* [(?*]] , (64) 
where we have used the expression (46) for m|- After implementing the RG-improvement, it becomes 

\ rg \r 

(M^o) 2 = M'=o + -§- [T[G*=o] ~ T[G*=o\] ~ [S[<2«=o] ~ S[G, =0 ]] - M| =0 . (65) 

It follows that, in the RG-improved case, the two definitions of the mass coincide at <j> = for any value of the 
temperature (as long as the masses are defined) whereas this was only true for T = in the non-improved case. 
The improvement has then restored a certain number of exact identities among the two possible definitions of the 



One can check that the corresponding differential equations are UV finite. This is not true if we would fix mo and A4 for instance. This 
is most certainly an artifact of the truncation. 
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mass. Similar remarks apply to the three different definitions of the four-point function at <p = and zero external 
momentum. In the RG-improved case they are identical for any temperature 

V$2o = V^ = ^=o. (66) 

The equality VTf? = V^—o is a particular case of the more general result 

V™ (K) = V^ ~ V^ Q [B[G^ ](K) - S[G 0=O ](O)] . (67) 
To obtain the latter, we start from Eq. (19) and apply the renormalization condition Vd,=o t* = A* to obtain 

V^o(K) = A* + A^ =0 (K) - A*(0) - ^ f G| =O (Q)A 0=O (Q) + ^ P G*(Q*)A*(Q*) . (68) 

We note next that, under the improvement, A^AT*) = A^o.t* as defined in Eq. (54) becomes equal to A^ (AT), 
from which it follows that the two integrals in (68) cancel after the improvement. Equation (66) is finally obtained 
by noticing that A* + A 0=O (A") - A*(0) = A* — \l[B[G ^{K) - £*[G*](0)]. Similarly, using Eq. (16) and the 
renormalization condition V^—o^ = A*, we obtain 

%=o = A* - \ f A 0=O (Q)G| =O (Q)^ =O (Q) + | A*(Q*)<^(Q*)K(Q*) ■ (69) 

From the definition of V*(AT*) in Eq. (55) and from Eq. (67), it is easily checked that, V<f, = o(K) and V*(K+) become 
equal under the improvement. Then, the two integrals in the previous equation cancel identically and V^f = V^—q. 

As we shall observe in Sec. IV, the critical exponents, which are of the mean-field type in the two-loop approximation, 
are modified after the improvement is considered. In particular, the exponent 5 gets closer (although it remains of 
the integer type) to its expected value in three dimensions. An unfortunate feature of the improvement is however 
that it is only defined in the symmetric phase: below a certain temperature T c < T*, which will be identified later 
with the critical temperature in the RG-improved case, the solution of the gap equation at vanishing field M^ = o is 
not defined. Therefore, it will be only possible to determine the improved critical exponents from above the critical 
temperature. In particular, we will not be able to access the improved value for the exponent (3. 



E. Multiply defined four-point functions 

In the next section, we solve the gap and field equations, using the expressions Eqs. (45), (46), (47), (50), (52) and 
(58) for the bare parameters, both in the two-loop and in the RG-improved two-loop approximations, and use the 
corresponding effective potentials to discuss the characteristic features of the phase transition in the model. Before 
we do so, however, it is interesting to study the temperature dependence of the three four-point functions V0=o, 
V<j,=a, and V^—o defined at vanishing field and zero momentum. Since G^ = o is a free type propagator, the four-point 
functions at zero field are all given in terms of perturbative sum-integrals. In this perturbative setting, it is then 
easy to check that the four-point functions are renormalized by the bare couplings Ao, A2 and A4 obtained in the 
previous section, without relying on the general proof given in Sec. VI. Moreover, since we are only interested here 
in the continuum limit, we can determine the renormalized four-point functions using any regularization. We shall 
use dimensional regularization or cutoff regularization, depending on our convenience. More precisely, V^, = q and 
V^—o, because they can be expressed solely in terms of tadpole, bubble and setting-sun sum- integrals (and their mass 
derivatives), see below, will be evaluated using dimensional regularization. In contrast, V0 = o cannot be completely 
reduced to these simple sum-integrals and we shall compute it using cutoff regularization. 

The renormalized expression for V^=o is trivially obtained by combining Eqs. (20) and (47) 

J- = f + \ [B[G*=o](0) - B*[G*](0)] . (70) 

Since 2?[G0 = o] is a one- loop integral involving a free-type propagator, it is clear that there is no divergence in this 
formula, as it can also be explicitly checked by a direct calculation, for instance using dimensional regularization. 
To obtain an useful expression for V$=o we consider its difference with V^=o- Using Eqs. (19), (20), and (39) and 
introducing the splitting A2 = A21 + <5A2 n i, we arrive at 



V^o - %= = A 2] + SX 2nl - A - XlB[G^o}(0) - %^ / G;L (Q) [a 21 + <5A 2nl - A - A^[G 0=O ](Q) 

1 Jq 
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1 H^ + aA2ni-A*B[G*=o](0) 



t4=o 



G$ =0 (Q) a 2nI -A^[G 0=o ](O) • (71) 



Using the expression for 6\2 n \ in Eq- (50) as well as the expression for A21 in the form of Eq. (59), we obtain 



V<t>=o 



6[G 0=O ](O) - B*[G*](0)]B*[G*](0) - ^ 



dS[G 0=o ] d5*[G*] 



dM; 



=0 



<im 2 



+ 



2A« 



(72) 



A direct calculation using dimensional regularization shows that the sum of the two square brackets in this formula is 
finite (see Eqs. (B12) and (B13)). As already mentioned, the four-point function V^—o contains a three-loop integral 
which cannot be reduced to simpler sum-integrals. We shall use the following expression: 



A 2 - 

t^,=o = A 4 + 3—=- (V^-o — An) + A + V0=o 

An 



A 2 I dS[G^ =0 ] 



2 ^ dMl =0 



A 2 dS[G^ ] 
An dMl 



4>=0 



K / G 2 =O (Q)^[G 0=O ](O), (73) 



obtained from Eq. (16) by using Eqs. (19) and (39). 
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FIG. 3: Temperature dependence of the three four-point functions V0=o, Vj, = o and V^,=o at parameters m,/T t 2 = 0.04 and 
A* = 3. In the case of V^,=o and V^=o, the lines are obtained using adaptive integration routines to evaluate their expressions 
derived using dimensional regularization. In the case of V$=o the line is obtained by evaluating the three-loop integral in the 
last term of Eqs. (73) and (60) as a double sum (N T = 2 10 and N s = 25 x 2 10 ), while in all the other integrals, including those 
in the bare couplings, the Matsubara sum are done exactly and the momentum integral are evaluated with adaptive routines 
at cutoff A/T* = 100. The points are obtained by evaluating the integrals as a double sums using N T = 2 x 2 10 non- negative 
Matsubara frequencies, N s = 13 x 2 10 values of the modulus of the the 3-momentum, while decreasing the cutoff A linearly 
from A/T* = 190 at T = T* to A/T* = 30 at T = T c . 




The variation with the temperature of the three four-point functions is presented in Fig. 3. Due to our choice 
of consistency conditions, the values of the four-point functions coincide at but in general they differ at other 
temperatures (with the exception of those values of T where two of the curves cross each other). As we shall see 
later, for those values of the parameters chosen here, the system undergoes a second order phase transition at some 
temperature T c . Above T Cl where the system is in the symmetric phase and the n-point functions are indeed defined 
at <f> = 0, we observe that V^—q and V^=o stay pretty close to each other which shows that the violation of the exact 
identity V^—q = V^—o is a mild one. The discrepancy is more important in the case of V^—a although the latter remains 
positive as long as T > T c . Note also that none of the four-point functions vanishes at T c . 

In the non-improved two-loop approximation, the curves below T c should not be taken too seriously because in the 
broken phase the n-point functions should be evaluated at the nontrivial minimum of the potential and not at <j> = 
(which is actually a maximum of the potential below T c ). The reason why we are able to follow the four-point functions 
at 4> = below T c is that the curvature of the potential M? is different from the gap mass M? n . There is then 
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a range of temperatures T c < T < T c where, although the curvature turns negative, the gap mass remains positive 
making then possible the evaluation of rt-point functions at = in the broken phase. This range of temperatures 
becomes more interesting in the RG-improved two-loop approximation because in this case the critical temperature is 
T c and not T c and it makes sense then to follow the four-point functions at = down to T c . In fact, since the three 
four-point functions become equal in this approximation, see Eq. (66), only the curve of V^—o is relevant. We note 
then that, in the RG-improved two- loop approximation, the four-point function vanishes at the critical temperature 
T c - As we shall see, this is directly connected to the modification of the exponent S. 

As a final remark, let us point out that studying the three four-point functions in the range T c < T < T c is interesting 
on numerical grounds for it gives valuable information concerning the discretization effects of the numerical method 
used to solve the model (see Sec. V for details concerning numerics). Namely, some integrals are infrared divergent at 
T c , and because of this V$=o goes to while V^=o and V^ =u diverges negatively. There is a competition between the 
square of the derivative of the setting-sun integral times V^q and the last term of Eq. (73), which both go as MZl 
as the mass M^ = o goes to zero, and this competition determines whether V^ = o diverges negatively or positively at T c . 
In order to obtain the correct divergence of V^—o numerically, a not too coarse discretization needs to be considered. 

IV. STUDY OF THE PHASE TRANSITION 

We now compute the effective potential and study how its shape changes as we lower the temperature T from the 
renormalization temperature T* down to T = 0. We shall first define the critical temperatures and evaluate them, 
then study the nature of the transition, followed by the thermodynamical observables and the critical exponents. All 
details regarding numerics are gathered in Sec. V, where we explain in particular how to accelerate the convergence 
of Matsubara sums and how to achieve accurate convolution routines. 



A. Critical temperatures 

We shall see below that the effective potential is convex at the initial temperature T*, with a single minimum at 
= 0. In other words, the unique solution of the field equation is = and the system is in the symmetric phase. 
Note that this result is not obvious a priori: at the temperature the renormalization and consistency conditions 
impose that the curvature of the potential is positive at = and thus that the potential is convex in the vicinity of 
= 0, but there is no obvious reason why the potential should be globally convex. 

As we decrease the temperature away from new extrema can appear, that is nontrivial solutions of the field 
equation. In particular, if a second order phase transition occurs at some critical temperature T c , nontrivial extrema 
are generated from = 0, because the curvature of the potential at = vanishes and turns negative. The critical 
temperature T c is then given by the equation 

M|=o,t c = 0. (74) 

This is an implicit equation for T c . Note however that since G^^o is a free-type propagator, the determination of 
T c only requires the calculation of perturbative sum-integrals. In this perturbativc context, it is also relatively easy 
to prove that the curvature at = and thus T c possess a continuum limit, without relying on the general proof 
of renormalizability that we give in Sec. VI. We start from the expression of the curvature at = obtained from 
Eq. (38) by using the expression (46) for ni2 ■ 

M|=o = m * + Y ~ T *[ G *]] ^ f [S[G^] - S^GA] . (75) 

Introducing the splitting A2 = A21 + 5\2n\ and adding and subtracting an appropriate term 15 we obtain 

Mj,„ = m? + ^[r[G„„l-r.[G,]]-^(Mj, -m;) 



15 The reason for adding and subtracting this term is that the sum of the second and third lines of Eq. (76) is finite, as it can be checked 
by a direct calculation or by using the results of App. A. 
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2 
G 



5[G^ = o]-5*[G*]-(Mi =0 -mi) 



2 ,dS*[G*] 



dm? 



(76) 



Now, using Eq. (59), we arrive at 



M. 



0=0 



^[r[G, =0 ]-r,[G,]] 



2 

A: 

2 

_± 
6 



M. 



= 



r[G 0=o ] -t*[g.] + (m 



0=0 



-^-y[T[G0 =o ]-7;[G*]] 
m2)iB*[G*](0)liS*[G*](0) 



A 



2 d<S*[G* 



5[G0 =o ]-^[G*]-(Mi =o -m^) 



dm 2 . 



(77) 



Using the gap equation (35) at <fi — and the expression (45) for m 2 ,, one sees that the first line is simply equal to 
M| =0 which is convergent. Moreover, a direct calculation using the perturbative expressions for the tadpole, bubble 
and setting-sun sum integrals, shows that the sum of the second and the third line is also convergent, see App. B. 
The convergent equation determining T c reads then 

A 2 r 1 
7^[G0 = o,t c ] — T*[G±] 



= M 



0=0, T c 



2 
6 



(M, 



0=0, T c 



m 2 )£*[G*](0) B*[G*](0) 



5c[G 0=o ,T c ] - 5*[G*] - (M| =0 , To - m*) 



dg*[G*] 
dm? 



with 



0=O,T C 



r c [G 0=o ,Tj-r,[G,] + (M 2 



=0,T C 



m 2 )2?*[G*](0) 



(78) 



(79) 



We solve this equation for T c for different values of the cutoff A. The continuum limit T c (oo) can be computed using 
any regularization. We find it convenient to determine it using dimensional regularization. 16 Both for T C (A) and 
T c (oo), the Matsubara sums are computed exactly and the remaining integrals over the modulus of the momentum 
are determined using adaptive numerical integration. These features allow for an accurate determination of T C (A) 
and also for the study of the convergence of T C (A) towards its continuum value T c (oo), as illustrated in Fig. 4. This 
represents a valuable information for the numerical resolution of the model since, due to memory limitations we 
cannot afford taking too large values of the cutoff while maintaining at the same time a good description of the 
infrared. We shall later use this accurate determination of T c in order to test our numerical code. 



Similarly, one can define a "critical" temperature T c for the gap mass, below which the gap equation at zero field 
has no solution [11]. It is defined by 

^0=o,f c =O. (80) 

The way T C (A) approaches its continuum limit T c (oo) is represented in Fig. 4. Note that when both T c and T c 
exist, one has necessarily T c > T c since the very existence of T c requires M|_ T to be well defined, according to 
Eqs. (78). However, these two temperatures do not coincide in general. 17 As already discussed in Sec. II B, one of the 
peculiarities of the 2PI formalism is that the different possible definitions of the two-point function do not necessarily 
coincide within a generic truncation. Because of this, in certain truncations, such as the two-loop approximation 
considered here, M| =0 is not equal to M| =0 and thus T c ^ T c . 18 In general, the temperature T c , where the curvature 
of the effective potential at </> = vanishes, either signals the vanishing of the field expectation value in a second order 
phase transition or it represents the lower spinodal temperature in a first order phase transition. The temperature 
T c resembles more the critical temperature in statistic physics, since the vanishing of the gap mass means enhanced 
fluctuations. Note also that below this temperature the potential at (f> = cannot be accessed because the gap 
equation does not admit a positive solution at = 0. 



In App. B, using dimensional regularization, we provide an explicit continuum version the curvature M^_ , see Eqs. (B10) and (Bll). 
The cases for which T c = T c correspond to m* = and A* > in which case T c = T c = T». If m* = and A* = 0, the temperatures T c 
and T c are not determined since the gap mass is identically zero for any temperature and the potential is identically flat. 
Note that there are other truncations such as the Hartree approximation, or the truncation that includes the "basketball" in &[<f>,G] 
which arc such that M? = M?_ f , and thus such that T c = T C . 
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FIG. 4: Cutoff dependence of the critical temperatures T c and T c (points) determined using numerical integration for parameters 
m*/T+ = 0.04 and A* = 3, and their convergence towards the continuum values T" c (oo) and T c (oo). The different convergence 
rate shown by the fitted functions (solid lines), 1/A for T c and 1/A 2 for T c , could be related to our choice of a sharp regulating 
function and the presence of a nonlocal sum-integral in the determination of T c . 




B. Nature of the transition 



The existence of a solution to Eq. (78) depends on the values of the parameters m%/T£ and A*. By determining 
those values of the parameters for which T c = 0, we can then separate the parameter space (mJ/T^A*) in two 
regions, corresponding to the white and grey areas depicted in Fig. 5. A point in parameter space for which the 
system undergoes a second order phase transition belongs necessarily to the white region, for which a value of T c 
can be denned. However, contrary to our discussion of the Hartree approximation [11], we cannot draw analytical 
conclusions on the nature of the transition in one or the other region. We thus select a certain number of points 
in each region and investigate the nature of the transition numerically. In Fig. 5 the points that we have explored 
numerically are indicated with a cross. 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 

m5/T5 



FIG. 5: The nature of the transition in a wide range of the parameter space inferred from the numerical study performed in 
the selected points indicated with a cross. In the two-loop approximation the boundary between the regions with no phase 
transition and second order phase transition is represented by the solid line. Along this line T c = 0. The dashed line is the 
boundary between the regions with no phase transition and first-order phase transition in the Hartree approximation. In the 
two-loop case this curve corresponds to T c = 0. The meaning of the other curves is the same as in Fig. 4 of [11]: the labels 
indicate the value of ln(A p /m*), where A p is the Landau pole. In the the region ln(A p /m*) > 5 our results can be considered 
cutoff insensitive for a cutoff A much larger than any physical scale but below the scale of the Landau pole. 
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We investigate first the nature of the phase transition for the parameters toJ/T^ = 0.04 and A* = 3. For these 
parameters the phase transition was of the first order in the Hartree approximation, as one can see in Fig. 4 of [11]. 
These values of the parameters were also used in [1] , where the model was solved in the same approximation, but in 
Minkowski space. Lowering the temperature from T*, we follow the value of the order parameter <p obtained from the 
solution of the coupled set of gap and field equations. We see in Fig. 6 that 4> remains equal to zero down to some 
value of the temperature, which turns out to be T c , since at this temperature the curvature of the potential at <f) = 
vanishes, as it is clear from Fig. 7. Below T c , 4> starts to develop a non- vanishing value. No other extrema appear 
between T+ and T c , as it can be cross-checked by looking at the shape of the effective potential. This is a radically 
different behavior than the one obtained in the Hartree approximation, where some nontrivial extrema appeared and 
became the absolute ones before the curvature at <fi = could vanish. 




T/T* T/T* 

FIG. 6: Left panel: The variation of the order parameter (f> with the temperature. Right panel: The temperature dependence 
of the first bin of the self-energy. The minimum of this curve coincides with the temperature value T c , where <j> starts to develop 
a nonvanishing value. The insets show the discretization effects near T c . The parameters are m* /T* = 0.04 and A* = 3. 



In conclusion, the temperature variation of the field expectation value </> shown in Fig. 6 corroborated by the change 
of shape of the effective potential with the temperature shown in Fig. 7 indicates a second order phase transition 
for the value of the parameters studied, in accordance with the result of [1]. We observe a similar behavior for two 
other points chosen in the white region and we believe that this result generalizes to a large part of the white region. 
It would be interesting to study our approximation for very small couplings and see if this region is dominated by 
the Hartree approximation, which would imply a first order phase transition, in line with the Monte Carlo results 
of [18]. This would be numerically challenging since for such small couplings, it would be difficult to differentiate a 
second order phase transition from a weakly first order one. For the point that we have tested in the grey region, the 
potential remains convex all the way down to T = 0. 

We can now go back to Fig. 5 and compare it to a similar figure obtained in the analytic investigation of the 
Hartree approximation in [11]. In this study the parameter space was also divided essentially in two regions. 
The separating boundary of the Hartree approximation is now represented in Fig. 5 by a dashed line and cor- 
responds to those points for which T c vanishes. In other words, for points below the dashed line there is no T c 
and the potential can be evaluated at <fr = down to T = 0. For points above the dashed line, there is a T c 
and the potential cannot be evaluated in the vicinity of </> = below T c . As the analytic investigation of [11] 
revealed, in the Hartree approximation, points of the parameter space above the dashed boundary correspond 
to systems which undergo a first order phase transition and points below the dashed boundary correspond to 
systems that remain in the symmetric phase. 19 The inclusion of the setting-sun diagram in the functional 
G] seems to change the nature of the transition to second order while enlarging the parameter space for which 
a transition occurs since the solid boundary line is pushed deeper in the no phase transition region of the Hartree case. 



The dashed boundary is in fact a very narrow band along the boundary line of the Hartree approximation and completely indistinguishable 
from this line at the scale of the figure. For points inside this band, although nontrivial extrema develop between T» and T = 0, the 
trivial minimum persists down to T = 0. 
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FIG. 7: The temperature evolution of the subtracted effective potential indicates a second order phase transition. The pa- 
rameters are m*/T* = 0.04 and A* = 3 and the discretization is characterized by A/T* = 100, N T = 512, and N„ = 3 x 2 10 . 



Similarly, we can study the effective potential using the RG-improved two-loop approximation. Note that since 
-^0=0 = ^4>-o (see Eq. (65)), the critical temperature which characterizes the sign change of the curvature at <fi = 
is also the critical temperature T c related to the vanishing of the gap mass. Since the gap mass at 4> = is not defined 
below T c , neither is the running of the mass in our improved scheme which is thus only defined in the symmetric 
phase. For those points in parameter space that we tested, we observed either no transition or a second order phase 
transition at T c in the sense that the potential remained convex down to T c where its curvature at </> = vanishes. 
In fact, the effective potential becomes very flat at T c . This is because, the fourth derivative of the RG-improved 
effective potential at <fi = equals V^— o,t and thus vanishes at T c . Of course we cannot test what happens in the 
broken phase since our RG-improvement is defined only for T >T C . Note also that the boundary between a second 
order phase transition and no phase transition coincides in this case with the boundary between a first order phase 
transition and no phase transition in the Hartree approximation. Note finally that the RG-improvement can also 
be applied to the Hartree approximation. In this case, one can show that the curvature of the potential for any 
value of the field changes from Mj = Mj + (V^ — \±)(j> 2 , where V$ is the generalization of V^—o to nonzero field, to 

(M£ G ) 2 = M 2 + (V,/, - V0=o)0 2 - It can also be shown that V$ increases with 4> 2 from which it follows that (M£ G ) 2 > 

for any value of </> and down to the temperature T c at which (M%2 ) 2 = M|-o = 0. Thus, in the RG-improved Hartree 
approximation, the potential remains convex down to the transition temperature. 



C. Thermodynamical properties 

In this subsection wc study the bulk thermodynamic properties of the model based on the pressure and quantities 
derived from it, such as the interaction measure (trace anomaly), the heat capacity and the speed of sound. 
The pressure is obtained from the subtracted effective potential given in Eq. (61) as 

p(T) = lim [A 7 (0)| - A 7 (0)j . (81) 

Actually, we cannot evaluate Aj((f>) exactly at T = because we can take into account only a limited number of 
Matsubara modes. Therefore, we are constrained to approximate the value of A 7 (^) at T = 0. In order to do so, we 
first determine Aj(<p) for smaller and smaller values of T by progressively increasing iV r , then we fit a high order 
polynomial to the available data points and obtain an estimation of Aj(<f>) \ T _ by evaluating it at T — 0. The value is 
accepted if the scaled pressure p/p SB is a decreasing function of T for T — > 0, where p S n(T) = tt 2 T 4 /90 is the pressure 
of the masslcss boson gas. In the opposite case, we increase iV T , determine the subtracted potential for smaller T and 
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FIG. 8: Bulk thermodynamic quantities as a function of temperature for two different coupling values: the scaled pressure 
p/Psb, entropy density s/ssb and energy density e/esb (upper row), the square of the speed of sound and the heat capacity 
(lower left panel), the trace anomaly (e — 3p)/T 4 (lower right panel). The insets in the upper panels show in a log- linear plot 
the perturbative results for the relative quantities obtained using a high temperature expansion up to and including 0(A«) and 
0{\l /2 ) (see Eq. (82)). The upper and left axis of the plot in the lower left panel correspond to c^, while the lower and right 
axis correspond to the heat capacity C. The discretization parameters are those of Fig. 7 and ml/T? = 0.04. 



redo the fit. 20 

Having determined the pressure as a function of the temperature, the energy density is given by e = —p + Ts, 
where the entropy density is obtained using a numerical derivative as s — dp/dT. The heat capacity C = de/dT is 
obtained numerically as the second derivative of the pressure: C = Td 2 p/dT 2 . The square of the speed of sound 
c 2 = dp/de is determined from c 2 = s/C and the trace anomaly of the energy momentum tensor T^ v is obtained 
as A = T^/T 4 = (e- 3p)/T 4 or equivalently as A = Td(p/T 4 )/dT. All these quantities displayed in Fig. 8 show 
nicely the second order nature of the transition: the scaled energy and entropy densities e/£ SB and s/s SB , and the 
trace anomaly display a cusp at T c , while the second derivative of the pressure with respect to the temperature is 
discontinuous, as displayed by the speed of sound and the heat capacity. The discontinuity is more pronounced at a 
larger value of the coupling. 

In the upper row of Fig. 8 we display the temperature dependence of the scaled pressure, entropy and energy 
densities calculated for two different couplings. These curves cross each other at the value of the temperature at 
which the scaled pressure has a maximum. This is because at this temperature, the interaction measure vanishes 
(since p/pss oc p/T 4 ) and thus e = 3p and Ts = Ap. Moreover, since — 3psB and Ts$ B — 4pg B , it follows that 
p/Psb = s/s SB = e/e SB . In the insets of these plots we compare the temperature dependence of the pressure with the 



For example at parameters mJ/T^ = 0.04 and A* = 3 we used A/T* = 100 and N s = 3 X 2 10 and increased N T from 512, used for 
T/T* > 0.2, to N T = 4 X 2 10 for 0.1 < < 0.2, and to N T = 8 X 2 10 for 0.06 < T/T* < 0.1. 
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first terms in the perturbativc expansion obtained at high temperature [34] 

Ppert(T) = Psb{T) 

where the neglected higher order terms depend on the chosen renormalization scale. Note also that, although the 
formula was obtained in [34] in the MS scheme, one can use it with coupling A* because the differences between 
the two renormalization schemes appear only at higher order in the coupling. The pressure obtained in the current 
approximation at coupling A + = 3 is closer to the 0(A*) pcrturbative result for T > 3.5T C . For the larger coupling 
constant, A* = 7, the pressure goes below the 0(A*) perturbativc result but for such high value of the coupling it 
makes less sense to compare to the perturbative expansion. 

At high temperatures the trace anomaly vanishes and e/(3p) goes to 1, in such a way that, interestingly, e — 3p 
is negative and its magnitude increases with the temperature. The fact that e/(3p) — > 1 is reflected in the square 
of the speed of sound, which approaches at high temperature the value 1/3, called the conformal limit because in a 
conformal invariant theory in three dimensions c 2 = 1/3. For low temperature the trace anomaly shows a bump, for 
both values of the coupling investigated. At the larger value of the coupling A* = 7 the cuspy structure becomes more 
prominent. These interesting features were already observed in [35]. 

D. Critical exponents 

There are six static critical exponents a, /3, 7, S, 77, and v, but, as a consequence of the static scaling hypothesis 
for the thermodynamic and correlation functions, which is verified in particular in the presence of a fixed point in 
the renormalization group flow [36], there exist four scaling relations between them, so that only two of them are 
independent. Usually 77 and v are chosen and the other exponents can be determined from 21 

a = 2-dv, /3 = (d - 2 + 77) ^ , 7 =(2-77K and S= ^ +2 ~ V . (83) 

2 a — 2 + 77 

Note however that there is a priori no reason why these relations should hold in a given approximation of the theory, 
such as for instance the two-loop ^-derivable approximation that we consider here. In what follows we determine the 
critical exponents in the two-loop and in the RG-improved two-loop ^-derivable approximations and discuss which of 
the scaling relations are fulfilled. 

1. Critical exponents in the two-loop approximation 

First of all note that there is a priori an ambiguity in the determination of certain critical exponents. For instance, in 
order to obtain the exponent 77, we should study the behavior of the propagator at criticality. One possibility is to study 
G. The corresponding critical temperature is T c and not T c and the propagator should be evaluated at 4> = down to 
T c . 22 But since G^, =u is local, we conclude that fj + = 0. We could instead consider the propagator obtained from the 
second derivative of the effective action, which generalizes the effective potential to non homogeneous configurations 
of the field. We would obtain a momentum dependent "curvature" 

M%=o{K) = KHZ + ml + ^ T[G 0=O ] - ^ S[G 0=O ] (K) , (84) 

where S[G^—o]{K) is the momentum dependent setting-sun sum-integral with propagator G<f> = o. At T c this self- 
energy is critical, in the sense that its value for K = vanishes. However, since G^ = n is massive, the corresponding 
propagator shows no anomalous dimension. We conclude then that t) + — 0. Then, even though the definition of 77 is 
ambiguous, in the present case, both approaches lead to the same result 77+ = i) + = 0, which coincides with that of 
the mean-field approximation. 



21 The first and third scaling relations are the Josephson and Fisher identities and instead of the second and fourth one one could use 
equivalently the Widom and Rushbrookc relations: 7 = /3(8 — 1) and a + 2/3 + 7 = 2 [37]. 

22 If one evaluates G at <j> = only down to T c > T c and at Gz for T <T C , G never reaches criticality, see Fig. 6. 
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Similar remarks apply to the critical exponent v. If we define the correlation length by £ oc M^_ , its scaling can 

be obtained by subtracting the renormalized gap equation at T c from the renormalized gap equation at temperature 
T, that is: 



m: 



A* 



=o 



[T[G*=o] ~ Tt o [G ] + M 2 =( A[G.](0)] 



with G (Q) = I/O 2 - 
justified since M^ = q 



Using the high temperature expansion of the tadpole sum-integral given in Eq. 
-> as T — !> T c , and neglecting the contributions of order M|_ , we obtain 
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(88), which is 
(86) 



from which it follows that u + = 1. We can alternatively define the correlation length from £ oc M^ Q . The way the 
curvature vanishes at T c is studied below when determining the exponent 7. We obtain that M^_ vanishes linearly 
as T — T c from which it follows that v + = 1/2. In order to solve this ambiguity, note that the nature of the transition 
is determined from the change of shape of the potential at T c . The relevant value for the critical exponent is thus 
v + = 1/2, which is again equal to the value obtained in the mean field approximation. 

The critical exponent f3 by fitting (f> to |T C — T\" . This requires first an accurate determination of T c from our 
numerical results. 23 We could proceed by locating precisely the temperature at which <f> starts developing a non-zero 
value. However, since the temperature derivative of cj> is infinite at T~ , it is easier to determine the value of T c by 
locating the minimum of the self-energy at the lowest available momentum and frequency: indeed the self-energy 
reaches a minimum value when <f> starts to develop a nonvanishing value. This is shown in the inset of the right panel 
of Fig. 6. Once T c has been determined the exponent /? can be fitted. As shown in the Fig. 9, the fit is compatible with 
the mean-field value /3 = 1/2. A similar method is used to determine the exponent 5. We introduce an external field 
h (this amounts to shifting the effective potential by —h<p), we set the temperature T to the numerically determined 
value of T c and fit <f> to h 1 / 8 . The results are compatible with the mean-field value 5 = 3, see Fig. 9. 
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FIG. 9: Numerical determination of the critical exponents f3 (left panel) and S (right panel) through fits to the field expectation 
value <j>. The parameters and discretization used are the same as in Fig. 7. The inserts show the dependence of the the fitted 
values of j3 and S on the upper boundary of the fitting window. As the fitting window shrinks, the fitted exponents converge 
to the mean field values 1/2 and 3, respectively. 



In order to obtain 7, we fit the susceptibility \ = d(j>/dh at h = to a power law \T C — T\ 1 . Because, in the exact 
theory 



dh 



h=0 



^7 



(87) 



This value of T c is not the same numerically than the one obtained from Eq. (78) using accurate numerical integration of perturbative 
integrals. We will later use these two different ways to obtain T c in order to test our numerical procedure. 
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we can also fit the inverse curvature of the potential to |T C — T| - ^. Note that in a given truncation, such as the 
approximation considered here, there is an ambiguity in the determination of 7 because there is no reason a priori 
why 7 should equal 7. Our numerical results for 7 are again compatible with the mean- field value 7_ =7+ = 1, 
see Fig. 10. Note that 7+ was obtained using dimensional regularization. Indeed, as we already discussed, in the 
symmetric phase, the formula for the curvature at (f> = involves only perturbative integrals which can be evaluated 
using dimensional regularization. Of course, since the curvature is finite, its continuum result does not depend on the 
regularization chosen to obtain it. The determination of 7 + and 7" would be numerically more involved. 
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0.998 
0.997 
0.996 
995 

1e-07 1e-06 1e-05 0.0001 0.001 0.01 0.1 
Upper boundary of fitting in |T - T c |/T* 

FIG. 10: Numerical determination of the critical exponent 7- (lower curve) in the broken symmetry phase and 7+ (upper 
curve) in the symmetric phase by fitting a(T — T c ) 7± to the curvature. In the symmetric phase the momentum independent 
curvature is determined using dimensional regularization. The critical exponent converges to the value 1 when shrinking the 
size of the fitting window around T c . The parameters are mJ/T^ = 0.04 and A* = 3. 




Finally, the heat capacity has already been determined in the previous section together with other thermodynamical 
quantities, see Fig. 8. It presents a discontinuity at T c , as it is the case in the mean-field approximation. To this 
behavior, one attributes conventionally the value a — for the critical exponent a. To summarize, in the two- 
loop <I>-dcrivablc approximation, the critical exponents coincide with those in the mean field approximation. In a 
sense, although it predicts the correct order of the transition, the two-loop approximation is not enough to produce 
non-analyticities in the effective potential which would modify the Ginzburg-Landau picture. 



2. RG-improved critical exponents 

As already explained, the RG improvement that we have introduced is only applicable for temperatures above the 
transition temperature which in this case is equal to T c . We can only determine the critical exponents by approaching 
the transition from the symmetric phase. In particular, we cannot access the exponent f3. 

An interesting feature of the RG-improved approximation is that, in the symmetric phase, there is no dif- 
ference between M^° and M^ = o. The determination of v£ G is then not ambiguous and coincides with that of 
v + in the previous section. Then z/+ G = 1 which differs from the mean field value 1 /2. The value of 77 remains equal to 0. 

In order to determine the exponents S RG and 7 RG , we can take advantage of some simplifications which occur in 
the RG-improved field equation at T c . Remember first that the RG-improved equation is obtained by applying the 
replacements m t — > M<^ = o and A* — > V^—o- Because V^—o goes to zero as T approaches the transition temperature T c , 
we will be able to neglect a certain number of contributions. Moreover, since M^ = o goes also to zero, we will be able 
to use high temperature expansions for some integrals calculated in dimensional regularization. We use, in particular, 
the expansion of the tadpole 
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from which wc obtain 



as well as [38] 
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In order to obtain the RG-improved gap and field equations, we can apply the above-mentioned replacements in 
Eqs. (35) and (37). The expressions for the bare couplings m|, A21 and JA211I become 
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(91) 

(92) 
(93) 



where we have used the expressions (45), (46), (50) and (59) for m 2 ,, m 2 ., <5A2ni, and A21. Using Eqs. (89) and (90), wc 
find the following behaviors for these parameters as we approach T c : 
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A similar analysis can be done for A4 which is expressed in terms of Ao and A21 in Eq. (58). The last integral of (58) 
involves a three-loop sum-integral which we do not compute and whose high temperature expansion is not known to 
us. Therefore, we evaluated this integral numerically at constant temperature and found that its value goes as MZf 
as the mass M^=o goes to zero. Since this integral is multiplied by V£ =0 (T), it gives no contribution as T — > T c . For 
the other integrals the HTE is known. Using Eqs. (89) and (90) we arrive finally at 
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Using these replacements, one obtains the following field equation (coupled to the gap equation) in the presence of 
the external field h 



h = 



Mi 



— m 

3 



3 ^ 



-0 2 + r f jG^j 



(96) 
(97) 



where Gq(Q) = l/Q 2 . In obtaining these equations we have also used the fact that since the setting-sun sum-integral in 
the field equation and the bubble sum-integral in the gap equation are multiplied by V^—o their contribution vanishes 
at T c . This can also be checked numerically. For instance, in Fig. 11, wc show the flattening of the momentum 
dependent gap mass M%_ (K) as we approach T c due to the fact that the nonlocal contribution to the gap equation 
vanishes in this limit. 



The round bracket in the field equation (96) is just M 2 =Q = = 0. Using Eq. (45) with replaced by T c and m* 

replaced by to express m 2 ,, and Eq. (47) at the reference temperature T+ to express Ao, one obtains the following 
renormalized equations: 



h= -6 Ml*. 
3 



T % [GfcfJ - T % [Go] + m r 840^(0) 



(98) 



As h — > 0, <f> —> and thus M 2 = — > 0, which justifies a high temperature expansion. Using the first terms in the 



expansion of the tadpole sum-integral given in Eq. 



the gap equation becomes quadratic: 
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FIG. 11: Flattening of the momentum dependent self-energy at constant cj>, as one approaches T c in the RG-improved case. The 
model parameters are m+/T+ = 0.04, A* = 3 and the discretization is characterized by A = 50, N T = 2 x 2 10 , N s — 26 x 2 10 . 



where terms of order Mj j, /T 2 and higher were neglected. At lowest order, one can neglect the terms of order Mj j, 

and obtain Mzf c ~ 2O7r0 2 /(3T c ). Plugging this result into the field equation in Eq. (98), one obtains the analytic 
value 6 RG = 5. We shall not present the numerical determination of S RG , for the simple reason that the HTE is very 
accurate in the region of h used for the numerical determination of <5 in the case without RG-improvement (see Fig. 9), 
and thus the solution of the gap equation in (98) is very accurately approximated by the solution of the quadratic 
equation (99). The value <5 RG = 5 can be understood as follows: assuming that the potential admits a Taylor expansion 
around = 0, the field equation at T = T c and for small h, reads 

¥ + O(f) , (100) 

4>=o,t c 

where we have used that the second and fourth derivatives of the potential at <f> = coincide with M|_ and V<f, = Q. 
Using that M^_ Q j, = and V^—q f c = 0, and assuming that the sixth derivative does not vanish, we obtain h cx 5 

and then i5 RG = 5. Thus, although the RG-improvement ensures that V0=o,? c = 0' which is a necessary condition for 
S to be larger than 3, the two-loop approximation is not sufficient to generate nonanalyticities in the field dependence 
of the potential which would yield a noninteger value for S. 

The value of j£ G can be determined analytically with a similar calculation, since it is given by the way the 
curvature at zero M? Q T behaves as we approach T c . We have already seen that M^ = o oc T — T c . It follows that 
7+, = 2. The numerical determination of 7^ is again simple and does not warrant a presentation. Similarly, one can 
determine analytically 7+, and one finds "f£ G = J RG . 

Concerning the heat capacity, this can be determined numerically down to T c through the formula 
C = —Td 2 ^f(<f>)/dT 2 , by applying to the effective potential the method described in Sec. HID. Around T c an 
analysis based on the high temperature expansion reveals that the heat capacity behaves as ~ a+ + b + (T — T c ) 
with the constant a + = 27r 2 T c /15 and 6+ = 347r 2 T' c 3 /135, independently of the value of the coupling. Unfortunately, 
we cannot conclude on the value of cv RG because we do not know whether there is a jump in the value of heat capacity 
at T c . The only thing we can state is that the heat capacity does not diverges as we approach T c + . 

Note finally that in the RG-improved case the last two of the four scaling relations (83) are fulfilled with d = 3. 
The two other scaling relations cannot be checked for we cannot access a or j3. 
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V. GENERALITIES CONCERNING THE NUMERICAL IMPLEMENTATION 

The resolution of the gap and field equations or the evaluation of the effective potential, together with the deter- 
mination of the bare parameters, involve various sum-integrals. Local sum-integrals of the form 



V[f}= f f(Q) 
Jo 



and nonlocal sum-integrals 



C[f,g](K) 



f(Q)g(K-Q) 



(101) 



(102) 



in the form of convolutions. We now explain how these various sum-integrals are discretized in view of their practical 
evaluation. We also present a method which allows to increase the rate of convergence of the discretized sum-integrals 
towards their exact result, leading to a sizable improvement in accuracy or computational speed. 



A. Discretization of the sum-integrals 



Let us start with the nonlocal sum-integrals because, as we will see, this puts some restrictions on the choice of the 
discretization. Sum-integrals of this type will be evaluated using Fast Fourier Transform algorithms. A convolution 
such as (102) can be written as 



where we have introduced the Fourier transform operator T and its inverse J 7-1 defined by 



7-[/](m,„,q) = J dr J d 3 x e^"^ -/(r, x) , 
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(104) 
(105) 
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with ft = 1/T. The functions /(iw n ,q) that we will have to deal with are invariant both under iuj n — > —iui n and 
under rotations (they only depend on the modulus q = |q| of the momentum). It follows that their Fourier transforms 
J-"[/](t, x) are invariant both under r — > f3 — t and under rotations (they only depend on the modulus x = \x\). 
Similar remarks apply to the inverse Fourier transforms. Using these properties, we arrive at 
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2 J dq sin(qx) ( Tqf(0, q) + 2T^2 cos(cj„t) qf(iu} n , q) 



where we need only the Matsubara frequencies Lo n = lixnT with n > 0. We can rewrite this as 
T[f}. =2n{F c ® F S )[U] and ^[f], = ® , 



(106) 
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where the notation /. means that the function / is multiplied by the modulus of its 3d argument, for instance 
f m (iu!„,q) = qf(iuj n , q), and we have introduced 



j.p/2 oo 

2 / dr cos( Wn r) /(r) and ^[f]^) = T/(0) + 2T V co S ( W „r)/(i W „) , 
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as well as 
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Note that if f(q) is the 3d Fourier transform of a rotational invariant function f(x), that is f(q) = J d 3 a;/(|a;|)e ?qx , 
then /. = 27rJ r 5 [/.], and in turn /. = ^F^lf.]- 

We have thus reduced the evaluation of the convolution to the evaluation of sine and cosine transforms whose 
discretized versions (DST and DCT) can be performed efficiently using one of the variants implemented in numerical 
libraries. 24 These variants differ in the type of boundary condition used when the original data is extended in view 
of performing on it the discrete fast Fourier transformation. As explained in Appendix E of [40], for the rotation 
invariant part, we use a discretization which avoids potential singularities as x — > and q — > 0, in that it does 
not store on the grid zero momentum and direct space values, and which matches the boundary conditions of the 
DST-II and DST-III formulas for the sine and inverse sine transforms, respectively (note however that the method 
that we put forward in the next section allows to reduce considerably the appearance of singularities in the UV). 
In momentum space, the highest stored momentum is the cutoff A and the grid is defined as kj, = (k + l)Afc, with 
k = . . . N s — 1 and Afc = A/N s the lattice spacing in momentum space, while in direct space, the grid is defined as 
x s = (s + ^)Ax, with s = . . . N s — 1 and Ax the direct space lattice spacing satisfying AxAk = n/N s . We retain 
N T — 1 positive Matsubara frequencies and the static mode u> n = 0, so that the available Matsubara frequencies 
are w„ = (27r//3)n = nAui, with n = . . . N T — 1. The corresponding temporal grid is defined as 7* = (t + \)At 
with t = . . . N T — 1 and At the temporal lattice spacing such that At Aw = n/N T . One can see that with this 
discretization, the discrete version of the cosine and inverse cosine transforms appearing in Eq. (109) are the DCT-II 
the DCT-III, respectively. 



In order to write the discretized version of the nonlocal sum-integral (103) in a compact way, we first introduce a 
shorthand notation for the sequence of discrete sine and cosine transforms which acts on an N T x N s array in which 
we store the values of the Matsubara frequencies and the modulus of the momenta. We define the following forward 
and backward discrete transforms 



T N ^ N3 [f(t,s)](n,k) = BCT-ll t [BST-II s [f(t,s)](t,k)](n,~k), 
-^ T W/M)KM) = DC3T-m B [DST-m i6 [/(n,fc)](n, «)](*, a), 



(111) 
(112) 



where n, t = . . . N T ~ 1 and k, s = . . . N s — 1, and the array /[t][s] is denoted for simplicity as f(t, s). The index 
indicates the part of the array on which the transformation acts. Note also that ^jv T ,iV s and n s are inverse to 
each other up to a numerial constant: T^ 1 N [.Fat^a^ [/]] = 4iV 1 -JV s /. This comes from the fact that DST-III and 
DCT-III are the inverses of DST-II and DST-II up to factors 2Af T and 27V S respectively. With the notation above it 
is easy to see, by using Eqs. (106) and (107), that the discretized version of the convolution reads 



CN T ,N B [f,g](n,k) = 



k + 1 
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[(p + l)/(m,p)] (t, s) ■ F^ >Na [(p + l)g(m,p)] (t, a) 



(n,k), (113) 



where n, m, t = . . . N T — 1, k,p, s = . . ,N a — 1, and the prefactor c = T 2 Ar(Ak) 3 / (8ir 3 ) contains the dimensionfull 
quantities arising from the discretization of the integrals. 

Next, we turn to the sum- integrals of the local type. Having stored on the grid N T — 1 positive frequencies and the 
static mode w„ = 0, this sum-integral will be approximated numerically as 



Vjv T ,iv,[/] 



n= _JV T +l 
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d 3 q 
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(114) 



where the notation [...]jv a refers to some quadrature rule, in practice we use the trapezoidal rule. After the exact 
evaluation of the angular integrals, one applies the extended trapezoidal rule [41] for the integral over q in the interval 
[0,A]. To obtain the formula, we include first the zero momentum (not contained by our momentum grid) in the 



We use the routines of the Fastest Fourier Transform in the West (FFTW) library [39], which contain a factor of 2 in the formulas of 
the DST and DCT transformations. This is the reason for separating factors of 2 in Eqs. (106) and (107). 
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sequence of points on the abscissa. Then, having N s + 1 equally spaced points, we apply the trapezoid rule on the N s 
intervals (0, Afc), . . . , ((N s — l)Afc, N s Ak) and obtain explicitly 
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B. Increasing the rate of convergence of sum-integrals 

In this section we take advantage of the fact that the asymptotic behavior of the propagator G(Q) is exactly given 
by 1/Q 2 in order to accelerate the convergence of the discretized sum-integrals towards their exact result. In order 
to illustrate the method we consider the tadpole sum-integral T[G] first. The most straightforward way to compute 
the latter would be as 

T[G] ~ Vnt,n s [G] . (116) 

The error that one should expect from this type of approximation is studied in detail in App. C. It is shown in 
particular that the error related to the finite number of Matsubara frequencies is directly connected with the rate at 
which the summand G(Q) = G(iu) n ,q) approaches zero at large n. Then, if in one way or another we are able to 
reorganize the evaluation of T[G] in terms of sum-integrals involving summands which decrease faster than G(Q) at 
large n, we will certainly reduce the error. Consider then the identity 

T[G] = [T[G] - T[G*]] + T[G*] = / 5G{Q) + T[G*] . (117) 

JQ 

The second term T[G*] involves the free-type propagator G* and the corresponding sum can be computed almost 
exactly, 25 see App. B. The first term involves a Matsubara sum whose summand SG(Q) = G(Q) — G*(Q) decreases 
faster than G(Q) at large n. Then, if we approximate the tadpole sum-integral T[G] by 

T[G}~V Nt , Ns {5G}+T[G*}, (118) 

we obtain an evaluation of T[G] which is more accurate than (116) for the same number of Matsubara frequencies. 

The same strategy can be applied to the bubble sum-integral. We can of course use the straightforward approxi- 
mation 

B[G](K) ~ C Nt , Ns [G,G](K) . (119) 
But we can instead reorganize the calculation of B[G](K) first, according to 
B[G](K) = [B[G}(K) - B[G4(K)} +B[G i ]{K) 

= [ [G(Q)G(K - Q) - G±(Q)G*(K - Q)] + B[G*](K) 
JQ 

= f G is {Q)5G{K-Q)+ f T SG(Q)G(K-Q) + B{G it ](K) 

JQ JQ 

[G*(Q) + G(Q)]SG(K - Q) + B[G*](K) , (120) 



where we have used G(Q) = G*(Q) + SG(Q) as well as the change of variables Q — >• K — Q. The benefit of the 
last expression is that it involves a contribution B\G+\{K) which can be determine almost exactly, see App. B, and 
a contribution in the form of a convolution with an integrand which decreases faster in the UV than the original 
integrand. Our final approximation for the bubble sum-integral is then 

B[G]{K) ~ C Nt , Ns [G* + G, 5G] + B[G+]{K) . (121) 



The Matsubara sum can be performed analytically and the momentum integral can be computed numerically using accurate adaptive 
integration routines. 
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This is a better approximation than (119) for the same number of Matsubara frequencies. 



Finally, consider the setting-sun sum-integral S[G] = J Q G{Q)B[G\(Q) . The straightforward approximation would 



be 



Instead, we write 



S[G] 



S[G]~V NrtNe [GC NT , Ne iG,G}] 



G(Q) [B[G](Q) - B[G*](Q)] + / G(Q)B[G*]{Q) 
Q J Q 

G(Q)[B[G]{Q) - B[G4(Q)] + f SG{Q)B[G ir ](Q) + S[G*} 
Q JQ 



(122) 



(123) 



The first term involves a summand which decreases faster than the original one G(Q)B[G](Q). Moreover, the inner sum 
(and the corresponding momentum integral) appears as a convolution and can thus be treated efficiently using Fast 
Fourier Transform algorithms. In the second term, the summand decreases again faster than the original summand 
and it contains a factor B[G+)(Q) which can be determined almost exactly. Finally the last term 6>[GJ can be 
determined almost exactly, see App. B. Our approximation for the setting-sun sum-integral will then be 



S[G]^V Nt , Ns [GC n ^ Ns [G, + G,SG}]+V Nt , Ns [SGB[G4] + S[G*] 



(124) 



C. Optimized equations and bare parameters 



In this section we gather the different equations that need to be solved and the expression of the subtracted effective 
potential and its curvature at <fi — which needs to be evaluated and, using the operators Vm t ,n 3 and Cn t ,n 3 , we 
put them in a form which is ready for numerical implementation. The first equation to be solved is the gap equation 
which reads 



M\K) = m \ + ^ + ^V NT , Ns [SG} + ^[T[G,}-r4G,}} 
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2 C Nt , Ns [G* + G,6G] (K) 



b 2 [B[G4(K)-B4G4(0)], 



(125) 



where we have used the expressions for <5A2 n i and m§ and it is clear from Eq. (125) that these are computed almost 
exactly. The other bare parameters relevant for the gap equation are computed as 
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Note that we could compute these bare parameters almost exactly as well. However, the proof of renormalization of 
Sec. VI reveals that these parameters absorb divergences in Vn ti n 3 [SG\. It is thus natural to compute Ao and the 
outer sum- integral of A21 with the same N T and N s . In contrast, the inner sum- integral of A21 is computed almost 
exactly for it has to do with the last term of Eq. (125) which is determined almost exactly. The gap equation can 
be coupled to the field equation in order to determine the extrema of the effective potential. In the presence of an 
external source h, the discrctized form of the field equation reads 



h = $ {ml + y ¥ + y V Nr , N . [SG] + ^ [T[GJ - %[G*}] 



~[ V ^,N 3 [GC N ^ Ns lG* + G,5G}] +V N ^ Ns [6GB[G4]] - y[<S[G*] -<S*[G*]] 
where A2 is evaluated as A2 = A21 + 2A*(l — A^/Aq) and the bare coupling A4 is computed as 



Ao I 



G 2 (Q+)[B4G4(Q+) - B4G*}(0)]' 



(128) 



(129) 
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Here, again, the appearance of the difference of bubbles is due to the almost exactly determined last term of Eq. (125) 
(see the steps leading from Eqs. (152) to (155)), therefore it is natural to perform these integrals almost exactly in 
the expression of A4, too. The discrctized effective potential reads 



A 7 (0) 



2% + 2^ 1 dqq 



Tlnll- e-< /T ) -T*ln(l 
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V Nt . Ns [SG} + [t[g*]-t*[g*\] 



(130) 



with e* = \J q 1 + m*. The derivative Sj/8^> of the effective potential in the formula above stands for an expression 
similar to the right hand side of Eq. (128), but with <j) replaced by cf>, while the coupling Ao and A4 are computed as 
in Eqs. (126) and (129), respectively. Finally, the discretized form of the curvature at vanishing field is given by 



M: 



A 



A- 



=0 



m 
G 



I + -4-VN r ,N.[*G*=o] + 4 [T[G*] - TAG,}] 



Viv T ,iv, [G^ = oCjv t ,j\t,[G* + G^=o,5G^=o]] + Vjv t ,jv. [SG^n 0[G*]] — -i [<S[G*] — <S*[G*]] , (131) 



where G 0=O (Q) = l/(Q 2 + M^ =0 ), with M^ =0 obtained from 



M| =0 



Y v N r ,N. [^G 0=o ] + ^ [T[G4 T4G,}] . 



(132) 



D. On the iterative solution of the equations 



The solution of the gap equation (125) at fixed value of the field, as well as the solution of the coupled set of gap and 
field equations (125) and (128) are obtained using iterations. We illustrate now this iterative procedure in the case of 
the solution of the field equation at vanishing external field h and comment on the convergence of the algorithm at 
different values of the coupling constant. 

First, we determine the bare couplings Ao, A21, and A4 given in Eqs. (126), (127), and (129), and evaluate those 
perturbative integrals in the gap and field equations which are defined at temperature T+. In case the of Ao, the 
convergence of the Matsubara sum is improved as described in Appendix C (see Eq. (C12)). The explicit expressions 
of the integrals evaluated using the adaptive numerical integration routines of the GNU Scientific Library (GSL) [42] 
are given in Appendix B. The quantities determined in this way will not change during the iterative process which, 
after the initialization of the propagator G with G*, consists of the following two steps: 

1. In Eq. (128) the double sums and the perturbative integrals defined at temperature T are computed with the 
actual propagator G(iu> n , k). If the sum of </>- independent terms in the curly brackets is negative, then <f> ^ is 
expressed by equating the expression between curly brackets with zero, if the sum is positive, then the trivial 
solution <j) — is considered. 

2. Using the value of (f> obtained in step 1., M 2 (iuj n , k) is determined from Eq. (125), by evaluating the double 
sums and the integrals with the actual propagator G, then the propagator is updated with 

G updatc ^n,k) = (w* +* 2 +M 2 (iw n ,fc)) _1 . (133) 

These two steps are iterated until the procedure converges to the desired accuracy. We monitored the value 
of the propagator at the lowest available frequency and momentum and stopped the iteration when both 
|GW(0, Ak) - G( i+1 )(0, Ak)\ /G (i+1 \0, Ak) < 10~ 8 and ^ - \ /fl^ 1 ) < 1(T 8 were satisfied. When this 

algorithm is used to determine 4>(T) by changing the temperature, then the iteration starts with the converged propa- 
gator obtained at the previous value of the temperature. In some cases, like the determination of the critical exponent 
(5, one has to work with nonvanishing external source h. In this case <p is obtained in step 1. by solving a cubic equation 
in the field. 

This simple iterative procedure fails to converge for large enough A* (~ 9), because the corrections are too large 
in each iteration. In this case somehow, the first iterative step bring us out of the attraction domain of the solution, 
if such a domain exists at all. This mostly happens when solving the gap-equation at fixed value of the field. When 
both the gap and field equations are solved iteratively, the procedure eventually converges for even larger couplings, 
but the number of iterations increases with the value of the coupling. In order to increase the speed of convergence or 
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to achieve convergence at all, one follows the procedure used in [29] and modifies the value of the updated propagator 
with a weighted average between the old value of the propagator and the iterated value. In this case, in the (i + l)th 
iteration the following updating method is used 

+ (1 -a)G {i) {iuj n ,k), (134) 

where a £ (0, 1] and M 2 {iui, fc)^ +1 ) is calculated with G^ . For large values of the coupling, A* > 9, one needs a < 1 
to achieve convergence at m 2 /T 2 =0.04 and = 1. For small couplings, a — 1 ensures the fastest convergence, and 
in an intermediate range, 7.5 < A* < 9, the use of a ^ 1 increases the speed of convergence. 

E. Cutoff convergence, discretization effects and the role of improvements 

From the proof of renormalizability given in Sec. VI we know that our results should becomes insensitive to the 
cutoff when the latter is taken large. However, since the proof is based on certain arguments that we can only 
verify using some perturbative estimates, it is interesting to check the cutoff insensitivity numerically, within a given 
accuracy. To do so, we have to pay particular attention to the discretization because we have to ensure that the latter 
does not distort the physics neither in the ultraviolet nor in the infrared regime. This means that as we increase the 
cutoff we need a good resolution in momentum space, that is small lattice spacing Afc, and also enough Matsubara 
modes taken into account. This represents a challenge for the judicious use of the available memory. Note that the 
implementation of the numerical improvements described above helps in this respect for the same accuracy can be 
achieved with less discretization points. For the same number of discretization points, the improved code is more 
computer time demanding because it involves the accurate numerical evaluation of perturbative integrals. However, 
in order to reach the same level of accuracy, the non-improved code needs to be run with a higher number of finer 
discretization, which requires an increased computer time as well. 

In what follows we shall illustrate on some physical quantities the effect of the discretization related to the use of 
the fast Fourier transforms to compute the convolution integrals, the extended trapezoidal rule for the momentum 
integrals and the finite number of Matsubara frequencies. After showing that these discretization effects are under 
control we will show also that the quantities of interest calculated with our best algorithm converge with increasing 
cutoff. 



1. Discretization errors due to the use of FFT 

The discretization errors related to the use of the fast Fourier transformation for the evaluation of the convolution 
integral can be easily illustrated with the help of the exact three dimensional convolution 

/(0 G ^^^i arCtail (^)' (135) 

where G(p) — l/(p 2 + M 2 ). This simple example is also relevant for our four dimensional study because it corresponds 
to the contribution of the static mode at finite temperature. Even though the integral in Eq. (135) is convergent, it 
will be interesting to calculate it using the same regulator as the one used in the four dimensional case. Using similar 
techniques as in Appendix B, we arrive at 

caigiw - l, Gfa)Gfa - »> - sk£ dtkG ^ rt>%'+M o ■ < i3e) 

|p-?|<A 

which can be evaluated accurately using adaptive integration routines. 

We can use the two results (135) and (136) to benchmark our method for evaluating convolution integrals and 
also to test how the continuum limit is approached. The different ways of computing the three dimensional bubble 
integral are plotted in Fig. 12. Note first that the bubble integral Ca[G] in the presence of a cut-off deviates from the 
continuum result already for values of the momentum much below the cutoff: at p = A/10 the deviation is already of 
5%. The result of a naive convolution on the momentum intervall [0, A] using discrete sine transforms stays close to 
Ca[G] up to p ~ A/2 (interestingly, it is closer to the continuum result for larger values of the momenta but this is a 
numerical artifact whose sign cannot be controlled in general). Instead, if we use the improved formula Eq. (121) (in 



+ k 2 + M 2 (iw n ,k) 



(i+l) 
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FIG. 12: Comparison of different methods of evaluating the perturbative bubble integral in 3d to the exact result. For more 
explanation see the text. The mass parameters are M 2 = 0.01 and m\ = 1 and the cutoff used was A = 500 in some arbitrary 
units. For the FFT we used N s = 15 x 2 10 modulus of the momenta. 



three dimensions), we can reproduce Ca[G](p) on the whole range of available momenta, up to p = A. This is related 
to the fact that, in the improved formula (121), one of the functions to be convolved decreases faster in the UV than 
in the original convolution: 5G(q) ~ q 4 instead of G(q) ~ l/q 2 . The overall picture remains the same when the cutoff 
is increased. 



2. Discretization errors due to the use of a finite number of Matsubara modes 

In Fig. 4, the temperature T c for which the curvature at <f> = vanishes, was determined for different values of the 
cutoff, by evaluating some perturbative integrals accurately (after the Matsubara sums were performed exactly, the 
remaining integrals were performed using adaptive integration routines). We can use these values as a benchmark to 
test the accuracy of T c obtained using the discretized version of its defining equation and to illustrate the effect of the 
improvements on the numerical procedure. Here we focus on the discretization effects related to the use of a finite 
number of Matsubara sums using thee different levels of improvements. 

The unimproved code avoids the use of any adaptive numerical integration and uses instead the most straightforward 
discretization for the quantities appearing in the expression (64) of the curvature. The convolution is evaluated using 
fast Fourier transforms with the formula (113) and all the momentum independent sum- integrals are approximated 
with a double sum: a sum over a finite number of Matsubara frequencies and a summation over a finite number 
of the modulus of the momentum, using the extended trapezoidal formula according to Eqs. (114) and (115). The 
momentum dependent bubble integral appearing in the setting-sun integral (122) and the expressions (127) and 
(129) for the bare couplings A21 and A4 are evaluated as a convolution, cf. Eq. (119). The difference in the partially 
improved code is that it uses an accelerated Matsubara sum in the tadpole integral and in the bubble integral with 
zero external momentum appearing both in the expression of the bare quantities and in that of the curvature itself 
(see Eqs. (Cll) and (C12)). The fully improved code uses the type of improvement done in the partially improved 
case cf. (C12) only in the sum-integral appearing in the expression (126) of the bare coupling Ao, but, as a major 
improvement, it uses the subtraction procedure described in Sec. V which involves perturbative integrals evaluated 
using adaptive routines and leads to modified formulas, with improved convergence for the tadpole, as well as bubble 
and setting sun integrals, in which the convolution is applied to functions which decrease much faster in the UV than 
the original propagators. 

The results for T c obtained within these three levels of discretization are shown in Fig. 13. In the plot on the left 
the result of the fully improved code (points) shows very good agreement with the accurate result of Fig. 4 (line) . As 
shown in the inset the test of the convergence of T c to the continuum result as the cut-off is increased required the 
increase of N s . The discrepancy between the points and the cutoff result T C (A) is mainly due to the evaluation of 
the convolution integral with Fourier techniques. Although barely visible in the inset, this discrepancy decreases with 
increasing values of the cutoff and N s . The scaling used in the left axis makes possible a direct comparison of this 
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FIG. 13: Illustration of the effects of the improvements by comparing T c determined numerically at parameters mJ/T* = 0.04 
and A* = 3 using fully improved, partially improved, and unimproved codes (see the text for explanations). Left panel: 
Solution of the defining equation (74) for T c , obtained with the fully improved code and with curvature M| =0 determined 
from the discretized expression (131) (points), in comparison with the solution of (78), obtained by evaluating the perturbative 
integrals with adaptive routines (solid line). Right panel: Dependence of T c on the number N T of Matsubara frequencies used 
to evaluate the sum-integrals in the expression (64) for M|_ in case of the completely unimproved code and of the partially 
improved code with A/T* = 100 and N s — 3 x 2 10 (inset). In case of the unimproved code, an asymptotic value of T c can be 
extracted with a fit, as shown by the line. 



figure to Fig. 2 of [1], where the same quantity was obtained by solving the model within the same 2PI approximation, 
but in Minkowski space. Note that in that reference T c (denoted there as T c ) slightly increases with the cutoff. This 
is not a shortcoming of the renormalization procedure, because here we have applied exactly the same method which 
leads to the same relations between the bare and renormalized quantities, rather it is probably a discretization artifact 
of the numerical method used in [1]. 

The plot on the right of Fig. 13 shows T c obtained with the unimproved and partially improved code (inset). The 
result of the partially improved code are acceptable if N T is increased by a factor of 5 compared to that obtained 
with the fully improved code. The values given by the unimproved code are far away from the true ones, even for 
huge values of N T . However, due to the decrease of the results with N T , an acceptable asymptotic value for T c can be 
extracted through a fit. 

The comparison presented in Fig. 13 shows that the acceleration of Matsubara sums discussed in Appendix C is an 
important ingredient of the numerical method used to obtain accurate results. 



3. Discretization errors due to the use of the trapezoidal rule 



The effect of the discretization related to the use of the trapezoidal rule to perform local type integrals can be 
easily seen by comparing the values of Ao and A2 evaluated accurately using adaptive numerical integration, with 
those obtained for a given discretization, that is for fixed values of N s and N T . The comparison is shown in Fig. 2. 
Since the acceleration of the Matsubara sum given in (C12) is implemented in the expression of Ao, N T does not play 
practically any role, and thus the comparison tell us up to which value of A the discretization in momentum space is 
acceptable. Based on this figure we concluded that N s = 3 x 2 10 is enough for A/T* = 100, but for A/T* = 500 it is 
not sufficient to obtain accurate results. 

A second example where one can see clearly the effect of the discretization of the momentum integrals is the 
variation of V<j, = Q and V^=o with the temperature from T* down to T c , as shown in Fig. 3. There we compared the 
values obtained using the numerical method used in the fully improved code with those obtained by evaluating the 
perturbative integrals adaptively. We saw that in order to be able to obtain for a given discretization the temperature 
dependence of V^—o and V$ = q with Fourier techniques, one has to decrease the value of A as one approaches T c , 
because as a rule of thumb a good description requires to have the lattice spacing in momentum space smaller than 
the propagator mass, that is Afc = A/N s < M^-o- 

As a third example, if one tries to determine T c from the discretized version of its defining equation in the fully 
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FIG. 14: T c obtained from solving (137) as a function of N 3 at two different values of A. The curve with the smaller A/N s 
converges faster. 



improved case 

ml + Y V *r,N, [SGo] + y [T[G*] - T,[G A] = 0, (137) 

where SGo = Go — G*, with the massless propagator Gq(Q) = 1/Q 2 , one runs into difficulties related to the fact that 
one cannot resolve the infrared behavior of the double-sum, which would require a momentum lattice spacing smaller 
than the mass scale. The best one can do here is to fix the value of the cutoff and increase N s , that is determine T c 
for smaller and smaller values of the lattice spacing in momentum space Ak = A/N s . The value of N T does not play 
a big role here, as we have tested by using N T = 2 10 and N T = 2 x 2 10 . As shown in Fig. 14, T c decreases as 1/N S . 
This allows to determine quite accurately through a fit the critical temperature, even from the discretized version of 
the defining equation. 

The variation with the temperature of the order parameter and of the first bin of the self-energy obtained with 
the fully improved code shows (see Fig. 6) that the discretization effects are under control in the fully improved code 
for fixed value of the cutoff. In Fig. 15 we show the cutoff dependence of the order parameter and the self-energy 
at different values of momenta at a given temperature. The quantities seems to converge as 1/A, and practically 
one could regard them as cutoff insensitive, to a good accuracy. As already mcntionned, a good description of the 
results at large cutoff values requires huge values of N s . Moreover, in order to see the scaling behavior with A, N T 
had to be increased for large values of the cutoff (A/T* > 200) as well. This is because the error made by cutting 
the Matsubara sums depends on A. In the case of the unsubtracted tadpole for instance, the error is ~ A 3 /N T T, see 
App. C. Increasing A without increasing N T would produce a cubic divergent which is not the correct UV behavior 
of the unsubtracted tadpole. 

All these tests convincingly show that the subtraction method described in details in Sec. V accelerates the Mat- 
subara sums and renders more efficient the evaluation of convolution using fast Fourier transformations. Therefore, 
it represents a reliable numerical method capable of providing accurate results. 

VI. PROOF OF RENORMALIZ ABILITY 

In this section, we show that the gap and field equations, and the effective potential, arc rendered finite by the bare 
parameters given in Eqs. (45), (46), (47), (50), (52) and (58). This result is nontrivial because the bare parameters 
do not depend on T or (f>, whereas the gap and field equations, or the effective potential, do. The proof is also 
the opportunity to illustrate some useful techniques which allow to discuss the ultraviolet behavior of Matsubara 
sum-integrals. 
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FIG. 15: Cutoff dependence of the solution to the coupled set of gap and field equations (125) and (128) at T/T* = 0.8 (points), 
scaled by the corresponding asymptotic value at A — > oo obtained by fitting (j)^ + c/A and M^,(fe, u) + c/t, w /A (lines) to the 
corresponding set of points. The convergence of M 2 (k ,u) is slower for higher momenta (right panel). The parameters are 
m*/T+ = 0.04, A* = 3. The discretization is characterized by Ak = A/N s kept fixed at the value 10/2 10 and N T = 512, with 
the exceptions of the points at A/T* > 200, for which, in order to achieve accurate results, iV T was increased by a factor of 2. 



A. Renormalization of the gap equation 

We discuss the renormalization of the gap equation first. Using the expression (45) for ibq, the gap equation (35) 
becomes 

M 2 (K) = ml + ^ 2 + ^ [T[G] - T*[Gy] - ^ 2 B[G](K) ■ (138) 

To proceed, it is convenient to decompose the momentum dependent mass M 2 {K) into a local and a nonlocal part, 
that is M 2 {K) = M- 2 + M^(K) with 

M 2 Em! + ^ 2 + ^ [T[G] - %[G A] , (139) 

Ml x {K) = -^ 2 [B[G](K) - E*[G*](0)] , (140) 

where we have used the decomposition A 2 = A 2 i + <5A 2n i and the expression (50) for 8\2 n \- Let us now discuss the local 
and nonlocal parts separately and show that they are both convergent. Using the results of App. A, the difference of 
tadpole sum-integrals appearing in M 2 can be written as 

T[G] - TAG*} = [ * 5G(Q+) + f 5a(Q) , (141) 
Jq* Jq 

where 5G(Q+) = G{Q±) — G*(Q*). Here G(Q*) is the analytic continuation of the Matsubara propagator G(Q) to 
Matsubara frequencies at temperature T*, where G(Q) is originally defined for Matsubara frequencies at temperature 
T. The second integral is a Minkowski- type integral over Q = (qo,q) with 5a(Q) = p(qa, q)e(qo)(n\ qo \ — n* qo \) and 

p(qo,q) the spectral density which enters the spectral representation (A2) of G(Q). As explained in App. A, this 
formula is useful to discuss the ultraviolet behavior of T[G] — 7^[G*]. Indeed, if we write 8G(Q+) = —(M 2 (Q+) — 
ml)Gl(Q.) + <? r (Q*) with <5 r (<3*) = (M 2 (QA - ml) 2 Gl{Q±)G(Q*), we obtain 

T[G] - %[GA = - [ T \m 2 (Q.) - ml)Gl(Q,) + Gr(Q*) + f Sa(Q) , (142) 
Jq* Jq* Jq 

where the first term is the only one that can generate divergences in the gap equation. The second term is a sum- 
integral at temperature T* whose integrand G I (Q is ) decreases fast enough at large Q irl see App. A. Moreover, the 
third term is convergent due to the presence of 5a(Q), see App. A. Note that the quantity M 2 (Q±) which appears 
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in the decomposition (142) needs to be regarded as the analytic continuation of M 2 (Q) to Matsubara frequencies at 
temperature T*. Plugging this decomposition into the expression (139) for M 2 , we arrive at 



M 2 = m l-^{M?-ml) 



G r (Q*) 



Sa(Q) 



(143) 



where we have used the separation of M (Q*) into a local and a nonlocal part. The first line is very similar to what 
appears when one considers the Hartrce approximation and can be treated along the same lines: dividing the equation 
by Ao, gathering the contributions proportional to M 2 — m 2 , using Eq. (47) and multiplying back the equation by A*, 
we obtain 
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(144) 



The first line is finite for both integrals are convergent, but the integral in the second line is still divergent. In order to 
treat this last integral, we need to discuss the nonlocal part M^(Q) first and then its analytic continuation A/ 2 j(Q*) 
to Matsubara frequencies at temperature T+. According to Eq. (140), the nonlocal part M^(Q) involves a difference 
of bubble sum-integrals which is shown to be convergent in App. A. We also obtain a formula for the analytical 
continuation M^K*) which is needed to complete the discussion of Eq. (144): 

B[G]{Ki) - = / *G*(Q*)[G*(Q* + JT*)-G*(Q* + £*)] 



5G(Q+) [2G*(Q* + K*) + 5G(Q+ + K*)] 



2 6a(Q)G(Q + K*) , 



(145) 



where we have introduced a general subtraction point L± for later purpose. This formula is used in App. A to show 
not only that the analytically continued difference of bubble sum-integrals, and in turn M 2 [(A'*), converges but also 
that it grows logarithmically at large A*, this logarithmic behavior being completely accounted for by the first term 
of Eq. (145) with L± = 0. But it is precisely this contribution which generates the remaining divergence in Eq. (144). 
This divergence is then T- independent and proportional to 2 , as it should if it is to be absorbed by A21. In fact, after 
plugging Eq. (145) with = into Eq. (144), we obtain 



M?=ml + ^ 



G r (Q*)+ / Sa(Q) 



2 V 



2 r T, 



A * A3 



5G(R+) [2G*(iJ* + Q*) + 5G(R, + Q*)] + 2 / 5a(R)G(R + Q*) 



(146) 



where the only divergent contribution, that is the sum-integral in the second bracket, combines with A21/A0 to yield 
1, according to Eq. (52). This completes the proof of renormalization of the gap equation. 

B. Renormalization of the field equation 

In order to prove that the field equation is rcnormalizcd by the bare parameters 777,2, A2 and A4, we prove a stronger 
result first, namely that the first derivative of the effective potential 



(147) 
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can be put in an explicitly renormalized form. Using the expression (46) for m\, we obtain 
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In App. A, we show that Eq. (see (A25)) 



T. 



S[G]-<S*[G*] = 3 / <JG(Q*)B*[G*](Q*) + 3 / fo(Q) B? [G*](Q) + A f 5 , 



(148) 



(149) 



where 5G(Q+) and 5a(Q) have been defined in the previous section, B^[G*](Q) is the retarded bubble contribution, 
obtained from (Q*) after analytic continuation and AfiS is a convergent quantity. From this decomposition, 

used together with A2 = A21 + A2 n i and Eq. (141), it follows that 



57 A 2 A * a c , A 4 j,2 , 1 



*G(Q*)A*(Q*) 
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MO)A^(Q) 



(150) 



where A*(Q*) was defined in Eq. (54) and A^(Q) is the corresponding retarded contribution. It is natural to replace 
A*(Q*) by V*(Q*) because V iI (Q ir ) is renormalized, as it is clear from Eq. (55). Similarly, it is convenient to replace 
&*(Q) by V^(Q). From the definition of V^ = o in Eq. (18) we can write 



V*(Q*) - A*(Q*) = V?{Q) - A*(Q) 
Plugging these formulae into Eq. (150), we obtain 
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(151) 
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SG(R i ,)+ / <5cr(7?) 



(152) 



Using Eq. (141) in the last term of this equation, we recognize T[G] — 7^[G*] which we can rewrite using the gap 
equation (138), analytically continued from K to Q*, as 



^| I r 



5a(R)\ = ^[T[G]-T4G*]] =M 2 (Q*)-ml-¥-[M-\lB[G\(Q ir )}, (153) 



where we note that the right-hand-side does not really depend on Q+ and can therefore be brought under the integral 
sign in the last line of Eq. (152). We obtain then 
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(154) 



where we have used that <5G(Q*) + (M 2 (Q*) - m^G^Q*) = G r (Q*). Using Eq. (145) with K* = L* = and 
Eq. (54) we arrive finally at 



$1 



a: 

6 



ml - ^ A f S 



G r (Q*)K(Q*) 



A 4 - 



A*(Q*)GJ(Q*)K(Q*) 



K(Q*)G*(Q*) 



5G(R+) [2G*(i?* + Q*) + 5G(i?* + Q*)] + 2 / 8a(R)G(R + Q„ 



(155) 



According to Eq. (53), the only integral which is still not convergent, that is the last sum integral in the first line of 
Eq. (155), combines with A4 to yield the renormalized result A*. This completes the proof of renormalization of the 
first derivative of the effective potential and in turn of the field equation. 
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C. Renormalization of the subtracted effective potential 

Let us first give a simple argument valid as long as T > T c , which includes the symmetric phase and part of the 
broken phase since T c < T c , as discussed in Sec. IV. In this range of temperatures the subtracted potential A 7 (</>) is 
defined down to d> = and we can thus write 



A 7 (0) - A 7 (0) 



(156) 



The second term in the right hand side of this equation is convergent from the discussion of the previous subsection. 
The first term is the contribution to A 7 (0) at vanishing field which coincides with that in the Hartree approximation 
and which is easily renormalizcd, sec for instance [11]. 

A more direct way, which is valid both in the symmetric and in the broken phase, 26 closely follows the derivation 
used for the field equation. From Eq. (153) we obtain 




M 2 (Q+) -ml-^-[X 2 - \lB[G](Q+)] 



Ml{Q)-ml-^[\ 2 -\lB R [G}{Q)] 



(157) 



where M R (Q) and Br[G](Q) are retarded functions obtained after analytical continuation. Using the equation above 
together with Eq. (145) with K+ = L* = and its analytical continuation in the subtracted effective potential (61) 
one obtains 



A 7 (0) 



7o (to*,A) - 7 £(m*,A) + - 



^ 4 J.2 



ln[l + E(Q)G*(Q)] -E(Q)G(Q) 
1 



1 



*G(Q*)A*(Q*) + - / MQ)Af(Q) 



2*66 



JG(Q*)E(Q*) 



5a{Q)t R {Q) 



8 



T, 



6G(Q+) [B[G](Q*)-S*[G*](Q*)] + / 5a(Q) [B R [G]{Q) - [G*](Q)] 



(158) 



where we introduced the shorthand notation E(Q) = M 2 {Q) — m 2 . In the square brackets the first integral is finite 
because the difference of bubble sum- integrals decreases at least as 1/Q*, while the second integral appears in the 
finite part of Eq. (A25). The combination in the round bracket is finite because it appeared in Eq. (150), which was 
proven finite. It remains to prove that the sum of the first two integrals is finite. To this purpose, note that since 
the divergences of these two integrals are overall divergences (this has to do with the fact that E(Q) grows at most 
logarithmically at large Q), they do not depend on the temperature. Moreover, since 



In 1 



1 



1 



S(Q)G*(Q)] -S(Q)G(Q) 

i*<2(Q*)S(Q*; = -- 

the two divergences are opposite to each other and thus cancel in the sum 



t 2 (Q)Gl(Q) + 0(t 3 (Q)Gl(Q)) 



t 2 (Q*)Gl(Q+) + 0(t 3 (Q+)Gl(Q*)) , 



(159) 
(160) 



Of course the potential is defined only for those values of the fields where the gap equation admits a solution. 
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VII. CONCLUSIONS 

We have studied numerically the temperature phase transition of the real (p 4 model from the two-loop ^-derivable 
approximation. Our analysis reveals that the inclusion of the setting-sun diagram in the 2PI effective action turns the 
phase transition into a second order type, which is believed to be the true nature of the transition in the model. The 
correct description of the order of the phase transition, which is also reflected in the behavior of some thcrmodynamical 
quantities, like the heat capacity, speed of sound and trace anomaly, represents an improvement as compared to the 
Hartrce approximation, where the transition is known analytically to be of the first order type. With this investigation 
we confirm (with a higher accuracy) and complete former numerical results obtained in Minkowski space within the 
same approximation [f]. 

Leaving the framework of strict ^-derivable approximations, we have checked that the phase transition remains 
of the second-order type even if the contribution of the setting-sun diagram is included only at the level of the 
2PI effective action, whereas the gap equation is solved at a lower (Hartrce) truncation, considerably simplifying in 
this way the numerical solution of the model. We will report on the rcnormalization of this "hybrid" ^-derivable 
approximation and some of its applications in a forthcoming publication [43]. 

In the present two-loop ^-derivable approximation, using a combination of analytic and numerical methods, we 
have also determined the static critical exponents which turned out to be of the mean-field type. Implementing some 
ideas borrowed from the renormalization group approach, we found that if we let the mass and coupling run with the 
temperature the critical exponents depart from their mean-field values (also they remain of the integer or rational 
type). Some of them can be determined analytically. In particular, the critical exponent S which characterizes the 
"magnetization" on the critical isotherm is found to be equal to 5, closer to its expected value in three dimensions. 
Diagrams with higher number of loops have to be included at the level of the effective action in order to have wave 
function renormalization and in turn a nonzero anomalous dimension. It is an interesting question to know what type 
of resummation will eventually produce critical exponents which are not integer or rational numbers. The 2PI-1/7V 
provides one example of such a truncation, see [31, 32]. 

Finally, the present work provides a concrete illustration, at finite temperature, of the general approach to renor- 
malization in the 2PI formalism developed in [28]. Following this approach, wc obtained expressions for the bare 
parameters which, in the present approximation, are given in terms of perturbative sum-integrals. Since we were 
interested in the effective potential and in thermodynamical quantities, we solved the approximation in the imagi- 
nary time formalism, which avoids the discretisation of sharply peaked functions such as the spectral density. The 
shortcoming of this approach is that we needed to determine the (slowly convergent) Matsubara sums numerically. 
However, owing to the simple asymptotic behavior of the propagator, we could accelerate the convergence of the 
Matsusbara sums. The same property allowed us to increase the accuracy of the momentum integrals which were 
computed using fast Fourier transform algorithms. The tests performed on different physical quantities concerning 
various discretization effects confirmed the gain in accuracy. It remains to be seen to what extent the numerical meth- 
ods developed here can be applied to more complicated truncations in the 2PI formalism, like the next-to-lcading 
order 1/N truncation in the O(N) model, the difficulty being that beyond the present two- loop approximation the 
asymptotic properties of the propagator are changed in a nontrivial way. 
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Appendix A: Thermal expansions 

When discussing the renormalization of the gap and field equations, we have to deal with differences of sum-integrals 
such as the following difference of two tadpole sum-integrals: 



T[G] - %[GA - / G(Q) - ( * G*(Q*) , 



(Al) 
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where G(Q) is the Matsubara propagator at temperature T and G+(Q*) is the free-type propagator l/(Ql + m 2 ) 
that we introduced in Sec. IIIC. The notations T* and are used to emphasize that the subtracted contribution 
7i[G*] differs from T[G] not only by its propagator but also by the different temperature entering the Matsubara 
frequencies: Q* = (iw*, q) with w* = 2-!rnT ir . In what follows, we bring Eq. (Al) and similar differences for the bubble 
and setting-sun sum-integrals to a form which is convenient for discussing their ultraviolet behavior. To this purpose, 
we make extensive use of the "analytic" propagator 

g( Z ) ee / ^ / +oc d J»P^2A , (A2 ) 

J qo Qo - z J_ 00 2tt q - z 

where Z = (z, q) and z belongs to the complex plane minus some possible points and segments of the real axis where 
the spectral density p is non zero. When evaluated for Z = Q, that is for z = iu n with u n = 2irnT a Matsubara 
frequency at the same temperature T than the spectral density p, we obtain the Matsubara propagator G{Q) which 
appears for instance in the first sum-integral 7~[G]. But we can also consider a "hybrid" propagator G(Q*) by 
evaluating the analytic propagator for Z = Q*, that is for z = iw* with ui* = a Matsubara frequency at a 

temperature different from that of p, and use this hybrid propagator to compute a hybrid sum-integral such as 
7~*[G] ee Jq* G(Q+). This type of hybrid sum-integrals will be useful in what follows. 



1. Tadpole sum-integrals 

Using the spectral representation (A2) for the Matsubara propagator G(Q) in T[G], and performing the Matsubara 
sum, it is a simple exercise to arrive at 

T[G] =11 p(q ,q)n qo . (A3) 

Jqo Jq 

Writing n qa = n* + Sn qo , we obtain 



qo •> q Jqo Jq 



T[G}= / / p(q ,q)n*+ / / p{q , q)5n qo . (A4) 



By repeating the calculation that leads to Eq. (A3), it is easily checked that the hybrid sum-integral 7~*[G] is nothing 
but the first term of Eq. (A4). Then, if we introduce the notations Q ee (qo,q) and 5a(Q) ee p(qo,q)5n qo , Eq. (A4) 
can be written finally as 

T[G}-T4G] = [ 6a(Q) , (A5) 



where the right-hand-side is a Minkowski-type integral, in the sense that it involves an integral over a real frequency 
qo. We can use this formula to express the original difference of tadpole sum- integrals (Al) as 

T[G] - T4G,} = [T[G] - T,[G]\ + [%[G] - T*[G*]] 

5a(Q) + f * 5G{Q+) , (A6) 



Q JQ, 

where we have introduced 5G(Q+) ee G(Q*) — G*(Q*) which involves the hybrid propagator G(Q*). As we now 
show, this formula facilitates the discussion of the ultraviolet behavior. Notice first that Sn qo = s(qo)5n\ qo \ with 
5n\ qQ \ ee n| go | — n f go |> an d thus Sa(Q) = p(qo,q)£(qo)8n\ qo \, where e(qo) denotes the sign function. It follows that 
the Minkowski integral in Eq. (A6) converges. Indeed, the integral over qo is cut off by 5n\ qa \ both for positive and 
negative values of go- Moreover, although we cannot really prove this fact, the spectral density is expected to decrease 
fast enough at large q and fixed qo: a perturbative estimate in the present approximation shows that the spectral 
density decreases like 1/q 8 at large q and fixed qo. As for the second term in Eq. (A6), it involves a sum- integral at 
temperature T*. If we write 

6G(Q+) = -(M 2 (Q+) - ml)Gl(Q+) + G r (Q+) , (A7) 

with G r (Q*) = (M 2 (Q*) — ml) 2 Gl(Q±)G{Q±) and use the fact that, in the present approximation, M 2 (Q+) is expected 
to grow at most logarithmically at large Q*, we see that the second term of Eq. (A6) diverges logarithmically and 
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that the divergence is generated entirely by the first term of Eq. ( A7) . For the purpose of renormalization, it is then 
convenient to rewrite Eq. (A6) as 

T[G] - %[G*\ = - [* {M 2 (Q*) - ml)Gl(Q*) + [ * G r (Q*) + / Sa(Q) . (A8) 
Jq, Jq* Jq 

We mention that the quantity M 2 {Q+) needs to be regarded as the analytic continuation of M 2 (Q), which is initially 
defined for Matsubara frequencies at temperature T, to Matsubara frequencies at temperature T*, just as G(Q+) 
which enters 5G(Q+) or G r (Q+) is the analytic continuation of G(Q). 



2. Bubble sum-integrals 

A similar analysis can be done for the difference of bubble sum-integrals: 

B[G]{K) - B*[G*](L*) = / G(Q)G(Q + K) - f * G*(Q*)G*(Q* + £*) ■ (A9) 

Jq Jq* 

Introducing an additional momentum integral by means of a 5- function in order to symmetrize the role of each 
propagator, using the spectral representation (A2) and performing the Matsubara sum, we obtain 27 

B[G](K) = I I [ I (2^6^ (p + q + k) p(q , q)p(p ,p) 1 ± ±?fe 

9o + Pa + iu 

(2^) 3 ^(p + q + k) p(q , q)e(q )p( Po> p) ^ 2 "' go1 , (A10) 

where, in going from the first to the second line, we have used the symmetry between (po,p) and (qo,q) as well as 
the formula 1 + 2n qo = e(q )(l + 2n.u i). Note that, in contrast to the case of T[G], it is not straightforward here to 
compare B[G](K) to the hybrid sum-integral S*[G](A^) because the external frequency has changed. We compare it 
instead to B^[G]{K) where the superscript (0) refers to the zero temperature Euclidean-type integral 

B^[G](K) = J^G(Q)G(Q + K) = J ^G(Q)G(Q + K). (All) 

Since the frequencies are not constraintcd in the zero temperature integral, we can consider the same external frequency 
as in B[G](K). We can now check that the contribution involving "1" in the second line of Eq. (A10) is nothing but 
the one we would obtain by computing B^°'[G](K). As for the contribution involving Jii go i, one can use the spectral 
representation (A2) to perform the integral over p$. After performing the trivial integral over p, one obtains finally 

B[G]{K) = B^[G](K) + 2 f a(Q)G(Q + K) , (A12) 

Jq 

where we have used the notation Q that we introduced in the previous section as well as <j(Q) = p{qo, q)e(qo)n\ qo \ . 
Note that the propagator G{Q + K) needs to be understood as the analytic propagator defined in Eq. (A2). Back to 
our original calculation (A9), we write as before 

B[G](K) - £*[G*](A0 = [B[G]{K) - B*[G](#*)] + [B*[G](F*) - B*[G*](L*)] , (A13) 

where we have introduced an arbitrary external momentum H+ for later convenience. The first bracket can be treated 
using Eq. (A12) and the second bracket can be put in the form of a single sum-integral at temperature We finally 
arrive at 

B[G](K) - B*[G*](L*) = f T * [G(Q*)G(Q* + HJ - G±(Q*)G*(Q* + L+)] 

Jq* 



It is convenient to use the parity of the propagator to replace G(Q + K) by G(—Q — K) at the beginning of the calculation. 
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+ / G(Q)[G(Q + K) - G(Q + #*)] 

JQ 

+ 2 f a(Q)G(Q + K)-2 [ <j*{Q)G{Q + #*) , (A14) 

JQ JQ 

where cr*(Q) = p(qo, q)e(qo)n* qg ^. This formula is the generalization of Eq. (A6) to the case of a difference of bubble 
sum-integrals with non-zero external frequency and momentum. We will also need its analytic continuation from the 
Matsubara frequencies at temperature T in K to the Matsubara frequencies at temperature T* of some vector K*. 
This continuation can be done readily on Eq. (A14) by replacing K by K*. The decomposition (A14) and its analytic 
continuation from K to arc particularly useful to discuss the ultraviolet behavior. Using G(Q*) = G ir (Q i ,)+SG(Q i ,) 
and choosing H+ = L+, Eq. (A14) becomes 

B[G](K) - B4G*]{L+) = [ G(Q)[G{Q + K) - G(Q + L*)] 

JQ 

+ f * 6G(Q+) [2G4Q, + L^) + 8G{Q+ + L*)] 
JQ* 

+ 2 [ a(Q)G(Q + K) -2 [ a*(Q)G(Q + L*) . (A15) 

JQ JQ 

Similarly, making the choice K = K+ = we obtain 

B[G](if*) - B*[G*](L*) = / *G*(Q*)[G*(Q* + A'*)-G*(Q i + L*)] 

JQt, 

+ / * 6G(Q*) [2G*(Q* + K+) + 8G(Q* + K„)] 
JQ* 

+ 2 [ Sa(Q)G(Q + . (A16) 

JQ 

Using similar arguments as those presented in the previous section, it is easily checked that all the contributions in 
Eqs. (A15) and (A16) are convergent. In particular the Minkowski integrals are convergent due to the presence of the 
functions a(Q), cr*(Q) or 8a(Q). The sum-integrals at temperature T* which appear in the second lines of Eqs. (A15) 
and (A16) are convergent because they involve either 8G{Q i! )G i ,{Q i , + K+) or 8G(Q+)SG(Q± + K) which decrease fast 
enough at large Q*. As for the first lines of Eqs. (A15) and (A16) they differ only by the temperature at which the 
sum-integral is computed and by the propagator which is used (T = and G(Q) for Eq. (A15) and and G*(Q*) 
for (A16)). Each of the integrands appears as a difference of two terms which decrease both exactly as l/(Q 2 ) 2 in 
Eq. (A15) or as l/(Ql) 2 in (A16). It follows that the integrands decrease strictly faster than 1/(Q 2 ) 2 and 1/(Q 2 ) 2 , 
and the corresponding sum-integrals are convergent. We will also need to know the asymptotic behavior of Eq. (A16) 
as K± becomes large. The Minkowski integral is subleading for it behaves as 

J_S<x(Q)G(Q + K*) ~ i f.HQ) ■ (A17) 

Some analysis that we shall not reproduce here allows us to argue that the second line of Eq. (A16) behaves also like 
1/K 2 , up to logarithmic corrections, and the first one behaves logarithmically. Thus, the dominant contribution at 
large K* is encoded in the first line of Eq. (A16). 



3. Setting-sun sum-integrals 



Let us finally consider the difference of setting-sun sum-integrals at zero momentum 

S[G]-<S*[G*]= / / G(Q)G(R)G(R + Q)- [ * [ * G*(Q*)G*(i?*)G*(fl* + Q*) . (A18) 

JQ JR JQ* JR* 

Introducing an additional momentum integral by means of a ^-function in order to symmetrize the role of each 
propagator, using the spectral representation (A2) and performing the Matsubara sums, we obtain 

S[G]= [ [ [ fff (2^\p + c i + r)p(r ,r) P ( qa , q )p(p a ,p) H ^ + n ~^-r (Alg) 

Jr Jqo Jpo Jr Jq Jp r + <?0 + Po 



4G 



Note that there is no singularity as the denominator 7*0 + qo + Po approaches due to the particular combination of 
statistical factors in the numerator. For later convenience, we shall then replace r + q + p by r + q + p a + ia. 
This does not bring any imaginary part because the sign of a can be chosen arbitrarily. Generalizing the approach of 
[44], we now write n qo = n* o + Sn qo and — ri— qo = —n*_ qo + Sn qo , and obtain 

S[G]= III III (2^(p + q + r)p(r Q ,r)p(q ,q)p(p ,p) << ~ tZ^V^o 

ro + qo + Po + ia 



r J q J pa J T J q 



3 / / / I I l(2nrs( 3 Hp + q + r)p(r ,r)p( qo , q )p{ Po>P ) Sn l { fl ^ 

ra+qa+Po + ta 



ro J qo J Po J r J q J p 



3/ / / / / I (2nf5W(p + < i + r)p(r 0l r)p(q ,q)p(p ,p) , (A20) 

ro + qa + Po + 101 



ro J qo J PO J r J q J p 



where we have exploited the permutation symmetry between the (po,p), (qo^q) and (tq,t) pairs. Applying the same 
steps that lead to Eq. (A19) to the hybrid sum-integral <S*[G], it is easily checked that the first line of Eq. (A20) is 
nothing but iS*[G]. Comparing the second line to Eq. (A10), we see that it can be written as Jq 6cr(Q)B+[G](Q + ia) 

where Q + ia = (qo + ia, q) and B+ [G] (Q + ia) = B^[G] (Q) is the retarded (if we choose a > 0) contribution obtained 
after analytically continuing B*[G](Q*). Finally, in the third line of Eq. (A20), we can perform the integral over po 
using the spectral representation which leads to the retarded propagator G{R + Q + ia), as well as the integral over 
p. We obtain finally 

S[G] =5*[G] +3 I 5a(Q)B*[G](Q + ia)+3 f f 5a(Q)5a(R)G(R + Q + ia) . (A21) 
Jq Jq Jr 

Once again, we can use this formula to compute the original difference (A18) of setting-sun sum-integrals for vanishing 
external frequency and momentum. We write 

S[G] - 5*[G*] = [S[G] - S4G}] + [<S*[G] - S*[G*]] 

= 3 / 6a(Q) B*[G](Q + ia) + 3 / / Sa(Q)S<r(R) G(R + Q + ia) 
Jq Jq Jr. 

+ I I* [G{Q*)G{R*)G{R* + Q*)-G*{Q*)G*{R*)G*{R* + Q*)}. (A22) 
Jq* Jr* 



This formula is the generalization of Eq. (A6) to the case of the difference of two setting-sun sum-integrals for zero 
external frequency and momentum. To obtain a formula suited for the discussion of ultraviolet divergences, we use 
G(Q*) = G*(Q*) + 5G(Q+). We obtain 



S[G]-<S*[G*] = 3 I * «<5(Q*)B*[G*](Q0 + / * / * 6G(Q i ,)8G(R i ,) [3G*(.R* + Q*) + 6G(R* + Q*)] 
Jq* Jq* Jr* 

+ 3 I Sa(Q)B*[G](Q + ia) + 3 [ f 6cr(Q)Sa{R)G(R + Q + ia) . (A23) 
Jq Jq Jr. 

It is also convenient to use Eq. (A12) to rewrite 

£*[G](Q + ia) = B*[G*](Q + ia) + [B*[G](Q + ia) - B±[G±]{Q + ia)] 

= B*[G*](Q + ia) + I [G(R)G(R + Q + ia)-G i ,(R)G„(R + Q + ia)] 
Jr 

+ 2 I a±(R)G{R + Q + ia)-2 f a*(R)G±(R + Q + ia), (A24) 
Jr Jr 

where u*(Q) = p ir (qo,q)e(qo)n* qo \ and p*(qo,q) = (27r)e((7o)<5(<Zo — I 2 ~ m *) i s the spectral density corresponding to 
G*(Q*). We can write finally 

S[G]-S*[G*} = 3 I **G(Q*)B*[G*](Q*) + 3 / fo(Q)B,[GJ(Q + ia) 
Jq* Jq 

+ /"*/* SG{Q±)5G(R+) [3G*(.R* + Q*) + dG{R+ + Q*)] 
Jq* Jr* 
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+ 3 / 5a{Q) [ [G{R)G{R + Q + ia) - G+{R)G*(R + Q + ia)] 
Jq Jr 

+ 6 f 8a{Q) [ a±{R)G{R + Q + ia)-6 [ Scr{Q) [ a*(R)G±(R + Q + ia) 
Jq Jr Jq Jr 

+ 3 / / Sa(Q)Sa(R) G(R + Q + ia). (A25) 
Jq Jr 

A similar discussion as the one presented in the previous sections shows that only the first two terms are divergent. 
The sum of the remaining terms will be denoted Af<S. 



Appendix B: Perturbative sum-integrals 



In this section, we give the explicit expressions for the perturbative tadpole, bubble and setting-sun sum-integrals at 
temperature T involving a free-type propagator G(Q) = l/(Q 2 + M 2 ) of mass M, both in dimensional regularization 
in d = 4 — 2e dimensions and in the presence of a sharp 3d cutoff for each propagator. We also discuss the monotonous 
behavior of the perturbative bubble sum-integral in the presence of a sharp cutoff. 



1. Perturbative results using dimensional regularization 



From Eq. (A5), in which we set = and Sa(Q) = (27r)niq 1 <5 (g§ — q 2 — M 2 ), we obtain the decomposition 



-Mr 



%[G] = Te {0) [G] + T e (1) [G]. T e (0) [G] is the tadpole at zero temperature, which is easily computed to be 

M 2 r 



T e (0) [G] = 



16tt 2 



1 M 2 
- + In — 

e 





M 2 


V, M 2 \ 


1 








6 32tt 2 





Ti- 



1 



(Bl) 



where we introduced the standard notation = 47r/i 2 e 7E with 7 E standing for Euler's constant, while "7v [G] the 
finite temperature part of the tadpole, given by 



r e (1) [G] = M 



2 c 



q n e 



2fi 



2c 



(27r) d -! e q ( 47r )^r(d=i) J 



dqq" 



(B2) 



with e q = \J q 2 + M 2 . Taking a derivative with respect to M 2 in the previous expressions, we obtain a decomposition 
of the bubble sum-integral at zero external momentum B 6 [G](0) = BP[G\(0) + BP[G](0) with 

M 2 7T 2 " 1 



and 



Bi 0) [G}(0) 



£«[G](0) = -,i 2 



167T 2 



1 , M 2 
- - In — 
e fi z 



32tt 2 



jd-l 



q d n e 



hr 



2 M 2 



0(e 2 ) 



i ci— 4 ^ e <i 



(B3) 



(B4) 



where, in this last integral, we have used the fact that d/dM 2 can be replaced by d/dq 2 in the integrand and we 

have integrated by parts assuming d > 2. We have expanded the vacuum pieces T}°\g] and B^[G](0) up to and 

including order e for later convenience. The thermal parts 7^ [G] and B^^G] will be needed only to order e° and 
we can thus take the limit e — > in those contributions. 

We proceed similarly for the setting-sun sum-integral. From Eq. (A21), in which we set = and also e r = 
y^(k — q) 2 + M 2 , we obtain the following decomposition: 



,4e 



SAG] = n 

+ 3^ 4e 



jd-l 1 



{2ir) d - 1 J {2n) d - 1 Ae k e q e r e k + e q + e r 



3/z 



4 c 



Jd-l 



k n F 



(2ir) d 2e k 



jd-l; 



d d l q n ek n £q 



(2tt)«*-i J (27r)^i 2e k 2e q 



B<°> ( £k + ia, k) + Bf\-e k + ia, k) 
1 



(e r + ia) 2 - (e fc + e q ) 2 (e r + ia) 2 - {e k - e q ) 2 



(B5) 
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The terms of the sum above contain in order zero, one and two statistical factors with positive argument and will be 
denoted respectively as iSi°^[G], <Se^[G], and 6>i 2 ^[G]. Note that in the present pcrturbative calculation, the regulator 
a plays no role if we assume M 2 > 0. This is always true for [G] for it is the zero temperature limit of a diagram 
which does not depend on a. For si^G], the analytically continued bubble contribution is evaluated on the mass 

(2) 

shell and therefore does not generate any imaginary part. Similarly, the denominators in <S e [G] never vanish since 
the equation = e 2 - (e k ± e q ) 2 implies A(k 2 q 2 - (k • q) 2 ) + 4(fc 2 - kq + <j 2 )A/ 2 + 3AI 4 = 0, which has no solution if 
M 2 > because the three terms are positive and one of them is strictly positive. 
We can also write the zero temperature contribution iS e [G] in a covariant form as: 



5i°)[G] = M 



At 



' d d Q 
{2n) d 



0- d G{Q)G{K)G{K - Q) . (B6) 



This integral can be evaluated with the method given in Sec. 11.5 of Rcf. [3]. Wc only have to obtain the 0(e) 
contribution of the integral J defined in (11.53) of this reference, which using our relation between d and e is J = 
3/e + 3 + e (3 - 11tt 2 /6 + 2*i(2/3)) + C(e 2 ), with *i(x) = d 2 T(x)/dx 2 being the trigamma function. Expanding in 



series of e one obtains 
5i°)[Gl 



(16tt 2 ) 2 



3 3 / M 2 3\ / / M 2 3V 5 5 2 \ (2 



+ Q(e), (B7) 



This expression agrees with the form given in [45] , upon exploiting a relation between specific values of the trigamma 

X 

function and the Clausen function C^x) = — J dO ln(2 sin (6/2)), namely 



2 



Since Be (fco,k) is Lorentz covariant (in dimensional regularization) , Be°^(±£fc,k) does not depend on k and can be 

pulled out of the integral in <s| [G]. In fact this constant contribution can be computed analytically. We obtain 
finally 

^ 1, l G l-i(7- 1 "^ + 2 -^ + '") 7:, " |G1 - (B8) 

where, for later purpose, it is enough to expand up to order e° the prefactor of T e [G\. Finally, the contribution 

(2) 

Se [G] will only be needed in the limit e — > where it yields a finite result due to the presence of the two thermal 
factors. After integrating over the angles, we obtain 

327r 4 7 e k J H e q A(k 2 - kq + q 2 ) + 3M 2 



Using these results, we can now check that the combination 



C(M,m+) = T[G] - %[G,\ + (M 2 - m 2 )6,[G,](0) B*[G*](0) - - 



S[G]-S4G,}-(M 2 -ml)^l 



, (BIO) 



which appears in Eq. (78) leads to a finite expression from which the scale /i drops out. To see this, we use that 
the expression inside the first pair of brackets is finite and that in terms of the form finite x divergent we have to 
expand the finite part to 0(e a ), where the integer a is the power in the most divergent, 0(e~ a ) piece of the divergent 

part. Note however that wc do not need to consider the order e term originating from 23i^[G*](0) in this first bracket 

because it is identically canceled by the contribution originating from 7*, e [G] in the one thermal factor contribution 
to dS+^G^/dm 2 . In the limit e — > the expression reads 
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FIG. 16: Region of integration for the bubble integral. 



M 



16tH 



16tt 2 [m 2 

T (1) [G] - T, (1) [G,] + (M 2 - m^PlGJfO) B^[G^{Q) 



M 2 \n^ ■ 



+ (A/ 2 - m 2 .) 



V3 



si 1} [GJ(0) 



dSj 2) [G,} 
dm 2 



Similarly, for the combination appearing in Eq. (72) 

2?(M,m*) = [B*[G*](0) -B[G](0)1b*[G*](0) 



dS*[G*] dS[G] 



dm 2 



dM 2 



(Bll) 



(B12) 



the following explicitly finite expression can be obtained: 



V(M,m+) 



1 



+ 



1287T 4 

1 

16tt 2 



2 m* m* 
M M 

M 2 



1 



16tt 2 



T (1) [G] T* (1) [G,] 



A/ 2 



ln^-2 + -) + ^>[G.](0) 



M 1) [G.](0)-i3 ti ^[G](0) 



1 

3 

(i)r 



dSi 2) [G+] dSW{G] 



dM 2 



(B13) 



2. Perturbative results using a 3d cutoff 



The techniques of App. A can also be applied in the presence of a 3d cutoff. From Eq. (A3) with p(qo,q) = 
(2-K)e(qo)S(q 2 — q 2 — m 2 ), we obtain 



TIG] 



2tt 2 



A 2 l + 2n £ 

dqq 



l 

8^ 



AeA — M 2 arcsinh 



1 

2T 2 



dqq 2 



(B14) 



where e q = \J q 2 + M 2 . In the second line of Eq. (A10), we can perform the integral over p trivially, as well as the 
integral over po using Eq. (A2). After some obvious shifts of the integration variables, we obtain 



B[G](K) 



p(9o,g)e(<7o)(l + 2n\ qo \)G{iuj - q ,k - q) 

qo ^|q|<A 
|k-q|<A 

l + 2n E V s 
— G{ilu - £ q , k - q) + G(iw + e q , k - q) 

|q|<A Z£ q 1 
|k-q|<A 



(B15) 



where in the second line we have performed the frequency integral. The angular integrals are computed using the 
geometrical picture shown in Fig. 16. With a 3d spherical coordinate system we integrate over the common region of 
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two spheres of radius A whose centers are separated by a distance k = |k|. The contribution is clearly if k > 2A. 
There is a non zero contribution if A < k < 2 A but we will be mainly interested in the case k < A. The integral over 
one angle gives 2ir, while the remaining angle 9 goes from to 7r if < g = |q| < A — k and from to arccos a < it if 



A - k < q < A, with a = (k 2 
formula 



k 2 )/{2kq) determined by the intersection "point" of the two spheres. Using the 



/ d(cos6)G(iuj±s q ,k-q) = f 

J a J a 



d(cos 9) 



2kqa + uj 2 =F 2iuje„ 



1 , k 2 
In 

2kq k 2 — 2kq + uj 2 =p 2iuje q 



its limit for a 



j 2 + k 2 =F 2ioje q — 2kq cos ( 
T, as well as \nz + lnz = In \ z\ 2 , the final form of the bubble integral reads (if k < A) 



(B16) 



B[G](K) 



1 



I6n 2 k 



A-k ^ 

dqq — 



A 1 
dqq — 



2n £ (k 2 + 2kq + lu 2 ) 2 
In 



4uT£ 



2 =.2 



9 

2n F 



(k 2 - 2kq + lu 2 ) 2 + 4w 2 ff2 



In 



(A 2 - q 2 



\uj 2 e{ 



(B17) 



Setting uj = and taking the limit k — > 
frequency and momentum 

B[G](0) = 



/A-fc t 9 (fc 2 -2fcq + W 2 ) 2 + 4w 2 £ 2 

0, one easily obtains the expression of the bubble integral at vanishing 



1 

8^ 



-A. 



1 + 2n F 



A 



dq ■ 



l + 2n r 



(B18) 



which can also be obtained from (B14) by taking a derivative with respect to M 2 and using an integration by parts. 

For the setting-sun sum-integral we start from Eq. (A20) with T* = 0. In this limit n* becomes — 9(— po) and 6n po 
becomes e(po) n | PO | • The first line of Eq. (A20) gives just the zero temperature contribution [G]. After performing 
the frequency integrals and the trivial integral over r, we obtain 



1 fl(A-|k-q|) 
|k|<A i|q|<A 4efe£ g £ r £k + E q + E r 



(B19) 



where we have renamed p as k and introduced e r 
the method given in [33] and one obtains 



M 2 . The angular integrals can be performed with 



S^IG] = 



1 



32tt 4 



A k 
dk- 

£k 



A-k 



dq ± ln £ -±±^ 

S Q Sk + £p 



A-fc 



£n £k + £ 



£A 



(B20) 



where e± = y/(k± q) 2 + M 2 . The second and third lines of Eq. (A20) can be treated simultaneously. We note first 



that in the second line the factor — 9- 



l>0 



1 + 9p g + 9 ro can be replaced by — 1 + 29 r 



e(ro), owing to 



the symmetry under q$ -H> r$ of the integrand. Combining this contribution with the third line of Eq. (A20) and 
performing the frequency integrals, we arrive at 



S[G}-SW[G}=3 [ f 0(A-|k-q|) 

J|kl<A J lal<A 



|k|<A J|q|<A 

Doing the frequency integrals, one obtains finally 



n eA 1 + n e q ) 

2e k e q 



1 



1 



(£k +£q) 2 



(£k-e q ) 2 



(B21) 



S[G] -5 (0) [G] 



32ir 4 



dk k 



Sk 



A-fc 

dqq — 



n €q 4(fc 2 + kq + q 2 ) + MI 2 



A 



dqq 



- n E 



In 



4(k 2 ~ kq + q 2 ) + 3M 2 
4(e q e k ) 2 - (A 2 - q 2 - k 2 



M 2 ) 2 



A-k 



M 2 (A(k 2 -kq + q 2 ) + 3M 2 ) 



(B22) 



3. Monotonous behavior of the perturbative bubble sum-integral 



Let us prove that the perturbative bubble sum-integral at finite temperature 



B[G](K) = T 



A e{A 2 -q 2 ) 



9(A 2 - (q-k) 2 ) 



(2tt) 3 io 2 +q 2 + M 2 K - lj) 2 + (q - k) 2 + M 2 



(B23) 
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decreases as one increases uj or k 
dependence, we consider 



with (we set k = k/fc) 

X X {K) = TJ2 

a 

1 2 (K) = TJ2 



|k|. This is clear for the frequency dependence. To prove it for the momentum 



dk 



B[G]{K) = 21 1 {K) + 21 2 (K) 



(B24) 



d 3 q 8(A 2 - q 2 ) 



6(A 2 -(q-k) i 



(2tt) 3 co 2 +q 2 + M 2 ((cv n - tu) 2 + (q - k) 2 + M 2 )- 



(q-k)-k, 



d 3 q 9(A 2 - q 2 ) 



<5(A 2 -(q-k) 2 ) 



(2tt) 3 w 2 + q 2 + M 2 (w„ - lo) 2 + (q - k) 2 + M 2 



(q-k)-k 



(B25) 



and prove that both T\(K) and T 2 (K) are negative. Let us treat the contribution T\(K) first. It is convenient to 
decompose the integration domain, which is defined by the O functions, into three different regions C, D and D, see 
Fig. 16. The region D corresponds to {q < A} n {|q — k| < A} n {(q — k) ■ k > 0}. The region D is the mirror 
symmetric of D with respect to the axis (q — k) • k = 0. The region C is {q < A} n {|q — k| < A}\(D U D). One has 
X\ = Ii + Ii + Ii . In region C (and also in region D), one has (q — k) • k < 0, from which it follows that < 0. 
In order to treat the remaining contributions, for each point q in region D, we introduce its mirror symmetrized 
q' = q'(q). More precisely, we have the formula 



q'(q) = q-2((q-k)-k)k. 



(B26) 



From this formula or from the geometrical interpretation given in Fig. 16, it is easily checked that |q' — k| = |q — k|, 
(q' — k) • k = — (q — k) • k and |q'| > |q|. Moreover the previous transformation has the Jacobian Jy = Sij — 2kikj/k 2 
which is such that J 2 — 1 and thus |det J\ = 1. We can now use this transformation to express the contribution from 
region D as an integral over region D. It follows that 



d 3 q 



(q - k) • k 



(2tt) 3 (K - oj) 2 + (q - k) 2 + M 2 ) 2 



1 



1 



; 2 +g 2 + M 2 uj 2 +q> 2 (q) + M 2 



(B27) 



which is negative because on region D, |q| > |q'| and (q — k) • k > 0. 

In order to treat the contribution X 2 (K) we can use the following geometrical argument. Noting that X 2 (K) 
represents the contribution to the variation of B[G](K) which comes from the modification of the integration region 
as k increases, that the two spheres separate apart in this cases and that the integrand of B[G](K) is positive, it 
follows that T 2 (K) is negative. Alternatively, we can also perform the angular integral using the 8- function. We 
obtain 



1 2 (K) 



167T 2 fc 2 ^— ' (w n 



A 2 + M 2 



dqq 



<1 



2 - k 2 - A 2 



A-k 



r 



M 2 ' 



(B28) 



which is indeed negative since q 2 — k 2 — A 2 increases on the interval [A — k, A] and is equal to — k 2 < for q = A. 



Appendix C: Rate of convergence of Matsubara sums 

In this section we study the rate of convergence of Matsubara sums and relate it to the asymptotic behavior of 
the summand at large Matsubara frequencies. Consider first the perturbative tadpole sum-integral T[G], which we 
approximate by 



v Nt [G] = t 

n~^N J ( 27r ) < + ? 



M 2 



In order to study how this finite sum converges to its limit T[G] 1 we introduce the error 

d 3 q 

\n\>N T ■ 



Z Nt [G] = T[G] -V Nt {G]=tJ2 J ^ 



) 3 L0 2 n + q 2 + M 2 



(CI) 



(C2) 
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Note first that a very simple bound of the error is obtained by setting q 2 + M 2 to in the previous expression. We 
obtain |£jv t [G]| < Co<p (N T ) with c = A 3 /(24tt 4 T) and <p (N T ) = J2\ n \>N T V" 2 - From 

2 1 2 

E ^ - 2 E ry = tv7 ' (C3) 

|n|>JV T n>N T y ' 

we obtain an even simpler bound \£n t [G]| < 2cq/N t . A numerical investigation reveals that the bounds are saturated 
already for values of N T < A/T, which shows that the bounds provide a good description of the error in this range of 
N T and in turn that the convergence of the Matsubara sum is slow. 

The functions CQipo(N T ) or 2cq/N t are in fact the first terms of asymptotic expansions of the error at large N T in 
some appropriate asymptotic scales which we now discuss. 28 We start from the following identity: 

1 K h k (-'\\ K + 1 h K + 1 

a + b ^-^ a K+1 a + b a Ji+i 

k— 

Note that, when a and b are both positive, \rk\ < Min (b K+1 / 'a K+2 , b /a that is the rest of the geometric series 
is bounded both by the last term kept in the sum and by the first term neglected. Plugging this identity into Eq. (C2), 
we arrive at 

K 

£n t [G] = ]T c k [G] Vk (N T ) + R<£> [G] , (C5) 

A:=0 

with 

^ = {in 2 )k +2 T 2 k+ i £ d ^ 2 (i 2 +M 2 ) k and Vk {N T ) = ^ . (C6) 

It can be shown that he functions <fk{N T ) form an asymptotic scale at large N T , in the sense that, Lp k {N T ) / <pk' {N T ) — > 
if k' > k. Now, from the inequality obeyed by rj», it is easy to show that 

\Rk\ < c K+1 [G]tp K+ i{N T ) . (C7) 

In other words 

K 

£n t [G] = J2 c kiG}<Pk(N T ) + 0& K+1 (N T )), (C8) 

fc=0 

for any value of K > — 1, if we decide conventionally that the sum is empty when K = — 1. This shows precisely that 
the error admits an asymptotic expansion in the scale of functions (pk(N T ). The functions ipk(N T ) admit themselves 
an asymptotic expansion in a scale of inverse powers of N T from which we could also deduce an expansion of the 
error in such a simpler scale. The advantage of (C8) is that the corresponding series converges as K — >• oo for fixed 
N T . We get then a better estimate of the error than with an expansion in a scale of inverse powers of N T for which 
the corresponding series turns out to be divergent. 

Let us now see how to use this results to accelerate the Matsubara sums. Note first that the coefficient cq[G\ docs 
not depend on M. The simple bound in (C3) can be used to accelerate the convergence of the Matsubara sum by 
writing 

T[G] = V Nt [G] + £ Nt [G] = [v Nt [G] + coMNt)] + [£n t [G] - cwo(iV T )] , (C9) 

and computing the first bracket instead of Vn t [G]. The new error is bounded by ci[G]tpi(N T ) which is parametrically 
smaller than co<po(N T ). Using that ipk(N T ) < 2/N 2k+1 and taking the dominant contribution of Cfc[G], that is 
c fe [G] - 2A 2fe+3 /((2fc + 3)(47r 2 ) fe+2 T 2fc+1 ), we obtain 



ci[G]vi(JV T ) < 3 ( A 



c [G}ipo(N T ) ~ 20tt 2 \N t T 



(CIO) 



28 An asymptotic scale is any collection of functions (ipk{N T ))k such that each tpk(^-r) does not vanish above some value of N T and that 
for k < k' , ipy (N T )/ipk (N T ) — > as N T — > oo. This last condition is also written ip k /(N T ) = o(^(AT r )). 
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It is easy to see that in combination with the trapezoidal rule (115) the improvement amounts to using the following 
approximate evaluation of the tadpole sum-integral: 

V Nt,N s [GJ = Vjv t ,jv s [GJ + — , (Oil) 

where the notation VN T .N 3 [f] was introduced in Eq. (114). Similar considerations for the bubble sum-integral implies 
the following improved evaluation: 

v [r2l _ v 2 (AkrNl + 2Nl 

V N„N. [Ct J - V NrtNa [G J + N 3 TS , ((-12) 

where we used that the analogous simple bound here is J2\ n \>N l/™ 4 < 2/(3iV^). We can do even better by noticing 
that the leading part of the coefficient Ck[G] is of the order or A 2fc+3 and it is thus independent of M. We consider 
then 

V^ T [G] = V Nt [G] - V Nt [G*] + T[G+] . (C13) 

where G* is a free-type propagator of reference for which we assume that 7~[G+] is known exactly. The error becomes 
now 

K 

£'n t [G] = ^]( c fe[G] — cfc[Go]) (fik(N T ) + 0(tpK+i(N T )) . (C14) 

k=l 

Note first that, because co[G] does not depend on M, the leading contribution of the error (k=0) has dropped. 
Moreover the leading contribution to the term of order k, which was previously of the order of A 2fc+3 is now of the 
order A 2k+1 M 2 . 

Of course the previous discussion is not very usefull in the case where the mass M is momentum independent 
because our acceleration method requires that we are able to evaluate T[G+] very accurately which is the same as 
computing the perturbative T[G] very accurately. The approach becomes interesting when we evaluate the tadpole 
sum-integral in the presence of a propagator G(Q) = 1/(G/ 2 + M 2 (Q)) with a momentum dependent mass. If the 
latter remains positive and grows only logarithmically at large Q, which is precisely what happens in the two- loop 
^-derivable approximation, then it makes a difference to evaluate T[G] from (CI) which has an error bounded by 
|£/v T [G]| < co<po{N) < 2cq/N t , or from (C13) which yields an error bounded by 

fe[G]-^[ G .|l<T £ / A "A . (015) 

\n\>N T V ' " 

which is expected to be suppressed with respect to coipo(N T ) if M 2 (Q) grows only logarithmically at large Q, as we 
verified numerically. 
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